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ABSTRACT 


The  modeling  of  damped  signals  as  the  impulse  response  of  a  pole- zero  system 
is  considered  for  a  broad  range  of  pole-zero  modeling  algorithms.  The  goal  is  to 
obtain  the  best  possible  fit  between  the  model  impulse  response  and  the  modeled  sig¬ 
nal.  Prony’s  method,  the  least  squares  modified  Yule- Walker  equations  (LSM'S'WE). 
iterative  prefillering.  and  the  Akakie  maximum  likelihood  estimator  are  com.pared 
on  known  test  sequences  for  a  variety  of  model  degrading  situations  le.g..  additive 
noise i  to  develop  an  understanding  of  which  methods  are  most  suitable  for  mod¬ 
eling  real  world  signals.  A  correlation  domain  version  of  interative  prefiltering  is 
al.so  introduced.  The  most  robust  algorithms  are  determined  to  be  LSN^'WE  using 
singular  value  decomposition  and  iterative  prefiltering  (including  the  correlation  do¬ 
main  version).  Modeling  several  laboratory  generated  short  duration  acoustic  signals 
confirmed  the  robustness  of  LSMVWE  and  iterative  prefiltering.  It  is  shown  that 
correlation  domain  iterative  prefiltering  outperforms  standard  iterative  prefilterine 
when  large  model  orders  are  required  for  accurate  modeling.  Shank  s  method  was 
determined  to  be  the  most  effective  method  of  determining  the  zeros  of  a  pole-zero 
model  when  a  tim^  domain  match  is  required. 
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I.  INTRODUCTION 


A.  RATIONAL  MODELING  OF  TIME  SERIES  DATA 

The  need  to  find  a  compact  parametric  representation  for  time  series  data  arises 
in  many  fields  of  study.  In  electrical  engineering,  finding  such  representations,  or 
models,  appears  under  such  topics  as  optimum  control,  optimum  filtering,  system 
identification,  model  order  reduction,  waveform  encoding,  and  spectrum  estimation. 
Other  fields  refer  to  such  modeling  as  time  series  analysis  or  forecasting.  In  digital 
systems,  the  model  parameters  take  the  form  of  difference  equation  coefficients  or. 
alternatively,  transfer  function  coefficients.  The  Piost  general  form  of  a  different  e 
equation  is  one  that  uses  both  feed-forward  and  feedback  coefficients.  The  corre¬ 
sponding  transfer  function  is  a  ratio  of  two  polynomials  in  the  comple.x  variable  r 
and  is  refered  to  as  a  rational,  or  pole-zero,  model. 

Historically,  the  more  general  pole-zero  model  has  been  used  in  relatively  few 
applications  compared  to  all-pole  (feedback  only  i  and  -ro  (feed  forward  only: 
motiels.  In  some  cases,  an  all-pole  or  all-zero  model  is  ...  most  appropriate.  More 
often,  however,  the  all-f)ole  or  all-zert)  model  is  chosen  because  the  optima!  model 
estimation  procedures  are  lietter  understood  and  easier  to  implement  than  those  for 
prjle-zero  modeling.  This  is  particularly  true  for  the  all-f)ole  casv  which,  in  many 
situations,  can  be  determined  fiy  solving  a  set  of  linear  equations.  Two  recent  (level 
opnients  have  led  to  increa.sed  activity  in  applying  pole-zero  models:  ( 1 )  lechnologu  a! 
ad  vances  in  digital  hardware  have  dramatically  reduced  computation  costs  and,  ('d! 
a  greater  variety  of  efficient  techniques  for  estimating  pole-zero  model  parameters  is 
now  available. 
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Literally  hundreds  of  pole- zero  modeling  estimation  algorithms  and  applications 
have  been  published.  The  majority  of  these  are  built  around  probablistic  or  stoch.ts. 
tic  modeling  techniques.  This  is  because  stochcistic  modeling  is  usuall\'  the  most 
appropriate  to  forecasting  [Fief.  1, 2.  3l  and  spectrum  estimation  [Ref.  4.  (i,  vs  here 

\ei\  little  i.--  known  about  the  system  input  which  produced  the  time  series  beim; 
modeled.  .-X  deterministic  methodology  has  also  been  used  in  which  the  input  and 
outjiut  time  series  are  available  and  the  linear  system  which  'best'  (usually  in  a  leas’ 
squares  sense)  produces  this  cause  and  effect  is  determined.  This  approach  is  usually 
fotind  under  the  topics  system  identification  [Ref.  7,  i''.  9.  10;  and  waveform  encodiiiE 
Ref.  111.  .-Xiiother  large  body  of  literature  which  exists  in  parallel  to  stochastic  and 
deterministic  modeling  is  that  of  reduced-order  modehm:  Ref.  12  which  larcelv  deals 
w;’;:  'lie  s'.’^tem  control  applications  of  pole. zero  modeling. 

.Mthoiitth  stochastn  an<!  <letermimst ic  pole-zero  mou<'i;ni:  h.tve  t!;e  s,i;:;e  goal 
and  use  tiie  same  mathematical  techniques,  the  di'-tmc'ion  i>  sicmlicant  in  the  way 
Ua’a  treated  and  iri  the  perlormance  cnterion  [lo-t ’nated.  specific, dly.  a’ lonar;' 
terpiirements  and.  as'-iim’it  icin'-  aliout  the  prol)a'‘Mi:ty  <iensity  functions  of  ratuioii;  pro 
,  es-e-  in  ‘•tof  hast  ic  modeling  can  be  limit  mg  wlien  deaiing  wit  li  real  world  s'cstemv.  In 
I  ot.trast .  ijeterminist ic  t echiiiqiu's  make  no  speciln'  icssiimpt urns  about  t  l,e  time  series 
to  be  rmcdeled  except,  of  course,  that  the  form  of  the  model  chosen  is  a[>propriate. 

One  particular  deterministic  nuxieling  problem  that  ha.'-  received  little  detailed 
attention  is  tliat  of  fiiiding  a  [lole-zero  mod«>l  when  the  tmie  senes  being  motii-Ie'i 
is  a  transient,  'mqnilse  response'-hke  waveform  Kxamples  of  situations  where  sin  h 
mofiels  tiiav  be  useful  include  m  mtalal  or  shock  analysis  of  mei  hanica!  sv  stems  [Ref 
!3.  1 4  .  antenna  response  to  elei  tromagnetic  pulses  [Ref.  I'l  .  and  wavelet  estima 
tKiri  m  siesrriic  studies  [Ref.  lf»i.  Role  zero  models  make  sense  for  transient  wave 
forms  because  such  a  model  is  sirurtnral  for  many  transient  signals  1  hat  is.  impulse 

■} 


response  like  transients  are  usually  the  result  of  a  system  of  oscillators  which  have 
been  excited  for  a  very  short  time  period  relative  to  the  natural  frequencies  of  the 
oscillators.  Pole  pairs  in  pole- zero  models  correspondingly  represent  the  resonant  fre¬ 
quencies  of  a  digital  system.  System  zeros  allow  the  initial  conditions  or  phasing  of 
the  resonant  frequencies  to  be  modeled.  This  combination  of  poles  and  zeros  allows 
signals  to  be  modeled  based  on  time  domain  matching.  When  an  effective  time  do¬ 
main  match  is  achieved  many  concerns  about  assumptions  during  modeling  become 
moot;  an  effective  time  domain  match  haus.  by  definition,  effectively  characterized  the 
signal  in  question.  i’re\'ious  work  has  concentrated  on  finding  onl\  transient  model 
poles  [Ref.  17.  18.  19.  20!  or  concentrated  on  one  particular  technique  of  pole-zero 
impulse  response  matching  [Ref.  21’. 

B.  THESIS  OUTLINE 

I  his  thesis  provides  a  performance  conifiarison  of  several  fiole-zero  modeling 
prni'f'fiures  ap, plied  t(j  the  proijlem  of  model  «-stimation  from  inifuilse  response  data. 

1  ne  jirocedure'.  compared  are  chosen  to  form  a  cross  scctio:!  of  o[)timali!\  and  com;  i- 
tational  <  om[ilexit\' of  aN'ailable  tecliniques.  In  Chapiter  11.  the  riK.idelinc  proce'diires 
■-ehTteci  for  study  are  described  in  detail.  1  he  nmiaining  three  chapters  are  concerned 
with  cfirnparing  the  performance  of  these  niodeling  techniques  and  are  organized  as 
follows; 

1  To  study  the  sjiecific  modeling  properties  of  each  method,  test  impulse  response 
seijiieiices  are  mofieled  in  Chapter  Thre«'.  lest  sequences  are  ronstrurt»‘d  to 
simulate  the  degradaf lotis  likely  to  be  present  in  'real  world'  transient  signals 
(e  g.  noise,  unknown  model  order!. 

2.  I  sing  the  results  obtained  in  Chapter  111.  laboratory  generated  acoustic  tran¬ 
sient  data  is  modeled  in  Chapter  IV.  This  data  is  considered  to  determine  the 


3 


performance  of  various  techniques  when  modeling  the  highly  complex  data  char 
acteristic  of  real  world  sources  about  W'hich  very  little  is  usually  known. 

3.  Chapter  V  summarizes  the  main  conclusions  drawn  from  the  results  in  Chapter 
Ill  and  I\’.  Recommendations  for  further  study  are  also  presented. 
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II.  POLE-ZERO  MODELING 


A.  OVERVIEW— THE  VARIETY  OF  TECHNIQUES 

Choosing  a  pole-zero  modeling  technique  from  among  the  many  available  tech¬ 
niques  can  be  difficult.  The  complexity,  applicability,  and  demonstrated  effectiveness 
of  the  different  methods  are  not  always  well  documented.  .Additional  consideration  of 
the  many  refinements  that  often  evolve  as  a  technique  is  applied  to  different  probh-nis 
can  lead  to  a  perplexing  array  of  tools  with  which  to  attack  the  modeling  problem. 
For  the  particular  case  of  modeling  by  impul.se  response  matching,  little  work  has  ap¬ 
peared  which  sorts  out  the  strengths  and  weaknesses  of  available  modeling  terfniiq’ie> 
or.  in  fact,  indicates  which  methods  may  or  may  not  be  applicable, 

I  he  basic  scheme  for  fitting  a  pole-zero  system  impulse  response  to  a  civeri  data 
sequeiice.  j-(n  i.  is  illustrated  in  f  igure  2.1a.  1  his  is  sometimes  referred  to  ;i>  the  dtrfcf 
model.  ,\s  we  shall  see,  the  formulation  of  this  problem  leads  to  a  set  of  nonlinear 
eqiiation<-  which  recjuire  the  use  of  iterative  technicpies  to  soKe.  To  overcome  t  lie 
'  ornp,exit ies  inherent  in  solving  nonlinear  ecjuations.  the  impulse  respor''e  matching 
prc.blem  cat!  be  ref(;rmulated  a.s  shown  m  Figure  2.1b.  This  ma\  bt'  referred  to  as  the 
indirrrt  method.  Note  that  these  two  formulations  are  not  equivalent:  the  error  of 
t he  direct  method  of  Figure  2.  la  is  c^fn)  =  x(u)-  h{n)  while  the  error  for  the  indirect 
method  of  f  igure  2.1b  is  c  ,{ri )  =  b(n)  —  a{n}»  x{u).  where  /i(  ri )  is  the  impulse  response 
of  the  pole- zero  system  Vieing  found.  6{n)  is  the  corresponding  sequence  of  numerator 
(■((efficients  and  ofri)  is  the  corresponding  sequence  of  denominator  roefficients.  The 
solution  of  the  indirect  problem  is  considered  suboplimal  in  the 
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Z*!  ri  ) 


-i  A  ( ; ) 
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Figure  2.1:  Two  pole-zero  modeling  problem  formulations:  (a)  direct  for¬ 
mulation  and  (b)  indirect  formulation. 


6 


I 

I 


when  the  modeling  error  goes  to  zero,  the  effect  of  minimizing  c,(n)  is  not  the  same 
as  minimizing  ej(nJ  in  the  direct  method. 

One  way  to  organize  pole-zero  modeling  techniques  is  shown  in  Figure  2.2.  In 
deterministic  waveform  matching,  the  ideal  equations  relating  the  proposed  model  to 
the  available  input  and  output  data  sequences  are  formed.  The  solution  which  ‘best 
satisfies  these  equations  is  chosen.  In  this  context,  ‘best  usually  means  the  solution 
which  minimizes  the  sum  of  squares  of  the  equation  error.  These  are  essentially  the 
Pronv  type  methods  [Ref.  17,  22]  (when  the  input  is  an  impulse  response)  and  lea-st 
squares  system  identification  [Ref.  7.  8,  ‘23]  methods  (for  general  input  sequences). 
.A  number  of  iterative  techniques  for  waveform  matching  have  also  been  proposed. 
Waveform  fitting  error  [Ref.  ‘24,  ‘25j  and  inverse  filtering  error  [Ref.  26.  27[  are  the 
criteria  most  used. 

Linear  stochastic  pole-zero  modeling  techniques  rely  primarily  on  estimates  of 
second  order  statistics  (auto-  and  cross-correlations)  to  estimate  model  parameters. 
Spectra!  estimation  has  been  a  driving  force  for  these  methods  which  are  based  on 
solving  some  form  of  the  modified  Yule- Walker  equations  (see  ('hapt<*r  llli.  Other 
methods  which  utilize  reflection  coefficients  [Ref.  28;  and  higher  order  statistics  [Kef. 
2‘L  'Ml  have  also  appeared. 

The  Tnanmuni  Ukfhhood  technique  seeks  parameter  estimates  for  which  the  ob¬ 
served  flata  is  the  most  probable  in  the  sen.se  that  its  conditional  probability  density 
function  (likelihood  function)  is  maximized.  This  technique  is  considered  to  be  sta¬ 
tistically  optimum  but  is  quite  difficult  to  use  because  its  implementation  generally 
requires  the  minimizaticjn  of  a  highly  nonlinear  function  [Ref.  1,  31.  32,  33].  Four 
widely  applied  methods  from  each  of  the  categories  of  Figure  2.2  will  be  used  in 
this  thesis.  A  number  of  improvements  which  have  subsequently  been  suggested  for 
these  methods  will  also  be  considered.  Most  recent  work  in  pole-zero  modeling  and 
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spectrum  estimation  has  occured  at  the  internal  boundaries  of  Figure  2.2,  i.e.  equiv¬ 
alent  linear  techniques  are  sought  that  perform  as  well  cis  modeling  formulations  that 
require  solving  nonlinear  equations  [Ref.  34,  35.  28,  36.  37],  The  application  of 
these  newer  techniques  to  transient  modeling  will  not  be  considered  here  since  we 
expect  that  the  performance  of  the  methods  chosen  will  in  most  cases  bracket  the 
performance  of  these  newer  techniques  with  the  possible  sacrifice  of  computational 
efficiency.  The  four  pole- zero  modeling  techniques  chosen  will  be  described  in  the  re¬ 
mainder  of  this  chapter.  The  transient  modeling  performance  of  these  methods  that 
will  then  be  compared  in  Chapters  III  and  IV. 


Linear 


Iterat  I  ve 


Deterministic 

t  '  t 

! 

F.quation  Error 
.Methods 

Waveform  .Match mi: 
Inverse  Filtermi: 

.'^tiK  lia.'itir 

Correlation  Equation 
Error  Methods 

i 

.Maximum 
Likelihood  Method^ 

Figure  2.2:  One  way  of  organizing  the  various  pole-zeru  modeling  tech¬ 
niques. 


B.  THESIS  MODELING  TECHNIQUES 
1.  Prony’s  Method 

One  of  the  best  known  indirect  techniques  for  matching  a  waveform  to 
the  impulse  response  of  a  linear  time-invariant  system  is  Prony's  method  [Ref.  22]. 
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This  method  is  in  fact  a  special  case  of  least  squares  system  identification  in  which 
the  system  input  sequence  is  taken  to  be  a  unit  impulse  and  the  numerator  and 
denominator  coefficients  are  determined  separately. 

The  pole- zero  modeling  problem  is  formulated  as  follows.  The  time  domain 

difference  equation  for  a  general  feedback,  feed-forward  system  can  he  written 

p  Q 

-  j)  =  Y^b,u(n  -  i)  (2.1) 

j=0  t=0 

where  u(ri)  is  the  input  sequence  and  x(ti)  is  the  output  secjuence.  When  the  input 
sequence  i.s  taken  to  he  a  unit  impulse  (unit  samph'  function)  and  the  output  is  taken 
tc;  h#‘  the  corresponding  impulse  response.  (2.1)  can  he  e.vpressed  in  tlie  form  of  a 


matrix  equation. 


■  x'Ot 

0 

0 

f  1  ■ 

1  k  ! 

i  :  1 

X(  I  ) 

x!0) 

0 

Ul 

h  1 

Xi  .\  -  1  i 

j'.V-2i  • 

•  JI.V  -  P)  _ 

L  . 

0 

wiiere  ,\  |v  the  numfier  of  data  points  used  and.  without  loss  of  generalization,  do  is 
'Ct  equal  to  one. 

I'.quatioi!  f2  2'  can  he  solved  hy  partitioning. 


where  a  =  il  (ij  ■  af>'J .  b  =  U\\  />,  ■  and  and  Xp  are  the  corresponding 

lower  ami  upper  partitions  of  the  data  matrix  in  (2.2).  I  lie  upper  partition  consists 


of  the  first  Q  -t-  1  rows  of  the  data  matrix  and  the  lower  partitioti  is  composed  of  the 


remaining  rows. 

The  scjlution  can  then  be  obtained  by  first  solving  the  lower  partition. 


Xaa  =  0. 


(2.4) 
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for  a  and  then  finding  b  from  the  upper  partition. 


XBa  =  b.  (2.5) 

If  A"  =  P  +  Q  +  1  then  (2.2)  may  have  a  unique  solution  and  the  model  impulse 
response  will  exactly  match  T(n)  for  0<n<P+4^+l.  This  is  referred  to  as  the 
Fade  approximation  [Ref.  22|. 

In  most  circumstances,  however,  the  length  of  the  availalile  data  sequence 
far  exceeds  P  1 .  It  is  then  desirable  to  us»‘  all  available  data  in  settinc  up  i  2.2 ' 

Ibis  leads  to  an  overdeterrnined  set  of  linear  equations  for  the  lower  partition.  No 
exact  solution  to  (2.4)  usually  exists  in  this  case.  The  relat lori'-hifi  in  ''2,4'  becomes 

X  .<3  =  e.  .  2.»' 

wliere  e  1-  the  equation  error  that  will  be  present  The  solution  of  .  2.4  ■  and  'J.ti' 
requires  )  part  it  lonmc  of  X^  as  follows. 

.  X.,  ;  X',  .  2  7 

where  is  the  first  column  of  X.^.  If  the  remaining  matrix  X'.,  is  square  and  of  full 
rank,  then  the  solution  to  (2.4|  is  given  bv 

a'  =  (X[^  r'x,.,  I  2.''  I 

where 

1 

a=  ,  .  i2.'.<i 

a 

Otherwise  the  /rasf  .s^uare.1  error  solution  of  (2.(1),  which  muiirnii'es  the  squared  error 
e^e,  is  given  by 

a'  =  X7x^  (2.10) 


10 


where  is  the  Moore- Penrose  pseudoinverse  of  X'^.  When  X'^  is  of  full  (column) 
rank  P.  the  psuedoinverse  is  given  by 

a'  =  (X'JX'^)-'X:,x,4.  (2.11) 


(See  IRef.  4.  pp.  2S-33]  for  an  example  of  this  derivation).  Otherwise,  the  {)seudoin- 
verse  is  defined  bv 


urx,4v. 

rr. 


(2.12! 


where  rr,.  u,.  and  v,  are  defined  by  the  singular  value  decomposition  i  represen- 

I  at  ion  of  X'., . 

»• 


X  4  =  Y. 


^2.13' 


t=0 

The  parameters  rr,  are  the  singular  values  of  X'.,.  u,  aiui  v,  are  the  corresponduic 
left  and  right  singular  vectors,  and  U’  is  the  rank  of  X',,.  Sec  'Ref.  3^  C}.  j-j'  fi.r 
a  more  detailed  explanation  of  singular  value  decomposition.  Once  a  is  known,  the 
ufiper  [lartition  of  (2.2i  can  be  solver]  by  smifily  carrying  out  the  matrix  mult '.nlnat  ion 
'ie-.rrii>ed  le.  the  left  hand  sirh- of  (2..')). 

2.  Modified  \'ule- Walker  Equation  Methods 

.•\  well  knriwn  class  of  pole-zero  modeling  technujues  is  based  of  soKing 
some  f<..rm  of  the  Modified  'tule- Walker  Erpiations  (.M'l'WKl.  I'hese  erpiatioti'  can  br- 
dr-veloped  by  multiplying  (2.1 )  by  x(n  —  k)  and  taking  the  expectation  of  both  sides 
1  his  vields. 


r  Q 

^r3*r„(n  -  k)  =  ^f>,r„j(ri  -  i)  (2.14) 

*=0  1  =  0 

where  is  the  autocorrelation  sequence  of  the  system  output  and  is  the 

cros.scorrelation  of  the  system  input  and  output,  if  the  original  input  to  (2  1)  is 
assumed  to  be  a  unit  variance  white  noise  sequence  then  the  cross  correlation, 
is  given  by 


II 


k 

=  YLHmi^k] 

k 


Equation  (2.14)  can  be  then  be 


written  as 


or  in  matrix  form. 


2_  —  k)  —  ^  b,fi(i 


*r=0 


r 

'rrll? 

r.rl  -  1  i 
^r./O) 

r^,  1  ~  P 1 

i  ,. 

1 

1 

1 

1 

r , .  1  Q  i 

~  li 

U-/': 

••  ^.-pCJ~}-p, 

'  r  !  1 

i  j 

^rj-l  ■+•  1  1 

j  ^  1 

e,,,  l>^ 

^  p  ^  1  , 

^rPQ  ^  P  -  \  \  .  . 
'  ^tPQ-P^ 

1 

'  i 

j 

L  '■'/'  J 

!  1 
b,/i  I  -  ;  1  ' 

=  '  j 


■J.K.i 


A, P,„„v.  . 

lined  separately.  Taking  a  lower 


result.s  in  !||f> 


partition  of  the  ]a,sf  /’ 


matrix  equation. 


^■rjuat  u»ns  j j)  '  L\  ]  7  / 


If  /'+  1  ««o™rr.Ut,„„  |,p  ^ 


12 


will  again  desire  to  extend  R/i  by  letting  the  index  n  in  (2.16)  run  beyond  Q  +  P  I 
resulting  in  additional  equations.  This  leads  to  a  an  overdetermined  set  of  equations. 

R^a  =  e  (2.19) 

which  will  not  in  general  be  satisfied  with  zero  error.  .As  before,  application  of  the 
psuedoinverse  results  in  a  least  squares  solution  of  (2.19). 

To  understand  how  the  MA'W’E  methods  can  be  used  to  match  a  time 
series  to  the  impulse  response  of  a  pole- zero  .system  observe  that  12.16)  describes  the 
relationship  depicted  in  Figure  2.3.  This  operation  can  be  equivalently  expressed  a.s 

r^^lnj  =  fi(T> )  •  h( —ri }.  (2.20) 

If  th'“  signal  to  be  modeled  is  assumed  to  be  a  pole-zero  systi'm  impulse  response,  then 
for  the  purpose  of  implementing  (2.17).  the  signal  being  modeh'd  ran  be  substituted 
for  hint  in  (2.20). 


Figure  2.3:  The  system  relationship  described  by  (2.16). 

The  equations  of  the  upper  partition  of  (2.17)  (the  first  Q  equations)  are 
not  linear  in  the  6  coefficients  and  are  generally  not  solved  directly  from  (2.17).  Once 
the  denominator  coefficients  have  been  determined  from  (2.18)  or  (2.19),  any  of  a 
number  of  techniques  are  available  for  finding  the  desired  transfer  function  numerator 
coefficients.  One  method  already  discussed  is  to  set  up  and  solve  the  upper  partition  of 
F’rony's  method  in  (2. .5).  Three  other  techniques  are  spectral  factorization.  Durbin's 
method,  and  least  squares  identification,  each  of  which  are  described  below. 
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a.  Spectral  Factorization 

Equation  (2.16)  describes  the  time  domain  relationship 

r„(/)*a(/)  =  b{l)*h{-l)  (2.211 

where  a{l)  and  b{l)  represent  the  sequences  of  denominator  and  nvinierator  transfer 
function  coefficients,  respectively.  Taking  the  r-transform  of  (2.21  i  yield.s 

c).-4(r)  =  1.  i2.22’ 


.Making  the  substitution 


//(--')  = 


.'Mc-'i 


in  (2.22)  and  rearranging,  results  in 


.4(c'’ )5„(;l.-T- 1  ==  Hiz'ilh: 


- 1 


r2.2:{i 


lo  utilize  12.24!  to  find  the  polynomial  coefficients  of  lh:\.  whnh  are  the  eieinent- 
! 'f  tiie  .■-c’quenre  hill,  we  rnust  perform  spectral  factorizaiiun  i»f  the  ‘-e.juer.re  result¬ 
ing  from  'he  cruivolut ion  of  the  three  secpiences  ut/l,  at  ■  /i.  and  r,,ii  .\  lieiaiieii 
exjilaiiaf ion  of  this  technique  can  be  found  in  [Hef.  or  Hef.  -lO 

b.  Durbin’s  Method 

Durbin’s  method  [Hef.  4lj  makes  use  of  the  property  by  wiiirh  a 
firoress  containing  zeros  ran  be  represented  b\’  an  all-pole  system  if  enougli  jioles 
are  used  The  first  step  is  to  filter  the  sequence  to  be  modeled  through  the  inverse 
of  the  previously  determined  all-pole  filter  coefficients  as  shown  in  Figure  2  -1.  I  he 
resulting  residual  .sequence.  .s(n).  will  nominally  be  an  all-zero  sequence.  .\  large 
order  all-pole  model.  Ai^rj,{:).  can  then  be  fitted  to  ,s(n)  to  obtain  the  relationship 
illustrated  in  Figure  2.6.  If  the  model  order  for  .4iarjr(')  is  sufficiently  large,  all  of 
the  ‘information’  in  .sfn)  will  be  contained  in  the  coefficients  of  If  an  all 
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Figure  2.4;  Filtering  process  to  generate  the  all-zero  residual  sequence  for 
the  application  of  Durbin’s  method. 


Figure  2.5:  The  residual  sequence  approximated  by  a  large  all-pole  model. 


pole  model,  is  then  ronstrurted  for  setpienre  cf  ro>  j]in>  in 

the  relationship  obtained  is 

t 

.1 .  I  ' )  ^  1 —  ■  ’  ‘ 

1  herefore  the  transfer  funrtion  1 /.-liarj.i  -  i  of  Figure  li  b  is  repla'-ei;  by 


/fir  '  5:  - - li  lin 

'.'.incli  yieifi^  the  desired  mo\'ing  average  lalbzeroi  niudel, 

c.  Shank's  .Method 

Consider  the  system  shown  in  Figure  '2A)  where  the  impulse 

response  of  the  pr«‘viously  <leterniined  alFpole  portion  of  a  pole  zero  mode!  iusinc 
Frony's  method  for  example).  Shank's  method  [Kef.  I'J'  is  to  satisfy-  the  relationshifi 
of  figure  J.fi  in  the  least  srpiares  sense.  This  relationship  ran  be  described  b\-  the 
rnatri.v  equation 


1.*) 


■  h^{Q) 

fiAiQ-l)  • 

f>x(0) 

■  bo  - 

f>AiQ  +  1  ) 

hA(Q)  ■ 

hAil) 

6, 

.  1) 

/m(.V-2)  ■ 

1 

■  -  1  -  Q)  . 

.  . 

x{Q)  ■ 

+ 

eB(Q+  1) 

.  -  1)  . 

.  e  s  ( A  -  1 )  _ 

H,,,b  =  X  +  eg.  CJ.’JN) 

Equation  (2.28)  ran  be  solved  in  the  manner  of  (2.6'  with  H  4  analogous  to  X'^  and  x 
analagous  to  x^.  The  bi,  coefficients  can  therefore  f>e  found  by  using  the  psuedoinverse. 


f  H'  "  I 


h  4 !  r( 


ffU'l 


Xl  Ti  I 


Figure  2.6:  System  for  Sfiank’s  method  determination  of  transfer  function 
numerator  coefficients. 


3.  Iterative  Prefiltering 


.■\n  iterative  l<‘chni()ue  fur  solving  the  direct  modeling  problem  of  Figure  2  1 
railed  itrrativf  prr  fill  f  rut  fi  has  be<-n  proposed  in  [Kef.  21i  .\n  effective  application  of 
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the  technique  has  been  reported  in  [Ref.  43).  The  presentation  below  follows  that  of 
[Ref.  40], 

In  this  method,  the  direct  modeling  problem  error  (see  Figure  2.1a). 

ed(n)  =  T(n)  -  h(Ti)  (2.29) 

is  expressed  in  the  alternative  form 

ed(n)  =  x(n)  -  b(n)  *  hA(n) 

=  x(n}  *  hA{n}  *  aln)  -  b{n)  *  b^in)  (2.30) 

where  a{n)  and  6(n)  are  the  sequences  of  transfer  function  denominator  and  numer¬ 
ator  coefficients,  respectively,  and  is  the  impulse  response  of  the  .AR  (all-pole) 

portion  of  the  model,  i.e. 

h^{n)  <==>  — ^ — 

.4(r) 

By  then  making  the  equation  error  iterative  (superscripts  represent  the  index  of  iter¬ 
ation), 

e'/Un)  =  jfn)  •  n)  .  a’+'frO  -  h'*Un)  .  h'Ju).  (2,31  ) 

the  least  squares  error  solution  for  a’'*''(ri)  and  6‘*'(n)  at  each  iteration  can  be  cal¬ 
culated  using  parameter  estimates  from  the  previous  iteration  In  matrix  form  (2,31  i 
becomes 

h\{P} 

x.jp^ii  x;fii  b\(p-^i) 

.  x,(.v  -  1)  x;,(.v  -  1  -  P)  /,y.v  -  1) 


b-AiP-Qj 
h\{P  -  Q^\  ] 

h\(S  -\-Q) 


1 1> 


(P) 

1) 


.'+1/ 


,\’  -  1 


1 


^0 


.32) 
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where  ih(n)  =  i(n)  *  h\{n).  Note  that  (2.32)  can  be  solved  in  exactly  the  same 
manner  as  (2.6). 

a.  Correlation  Domain  Iterative  Prefiltering 

Many  pole-zero  modeling  algorithm's  which  were  originally  conceived 
bcised  of  the  time  domain  pole-zero  difference  equation,  (2.1).  have  been  reformulated 
in  the  correlation  domain.  Examples  of  this  include  the  correlation  domain  extension 
of  Prony's  method  resulting  in  the  modified  Yule- Walker  equation  methods,  and  an 
instrumental  variable  method  of  least  squares  system  identification  which  Soderstrom 
has  indicated  is  simply  a  correlation  domain  formulation  of  least  squares  system 
identification  [Ref.  44].  In  modeling  trials  conducted  for  this  thesis  both  of  these 
correlation  domain  methods  were  found  to  be  significant  improvements  over  their 
time  domain  counterparts. 

Iterative  prefiltering  can  also  be  extended  into  the  correlation  domain. 
To  see  how  this  is  done  first  note  that  the  direct  formulation  of  the  pole-zero  modeling 
problem  of  Figure  (2.1a)  may  be  reinterpreted  in  the  correlation  domain  by  employing 
the  relationship  of  (2.20)  and  Figure  2.3.  F'igure  2.7  illustrates  this  new  direct  pole- 
zero  modeling  interpretation. 

Proceeding  as  for  the  time  domain  iterative  prefiltering  above  we  write 
the  correlation  domain  error  equation. 

fd.corrin)  =  r„(n)  -  r(-n)  • /»(n) 

=  *  hy^[n)  •  a{n)  -  x{-Ti)  •  h{n)  •  h^{u)  (2.33) 

where,  as  before.  a(n)  and  b(n)  are  sequences  of  the  model  coefficients.  h,^{n)  is  the 
impulse  response  of  the  AR  (all-pole)  portion  of  the  estimated  model,  and 
is  found  using  xin)  as  the  desired  impulse  respon.se  .so  that  r„(n)  =  j-(n)  •  j(-ri). 
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Figure  2.7:  The  direct  pole-zero  modeling  problem  formulation  expressed 
in  the  correlation  domain- 

W  hen  the  error  i?  made  iterative,  (2.33)  becomes 

=  •'■(”)  •  j(-r!)  •  h\{n)  •  -  j(-r))  *  h''^Un}  *  h\in)  (2.34) 

where  the  superscripts  represent  the  index  of  iteration.  In  matrix  form  (2.34 )  becomes 

VKiP)  ■■■  r;(0)  T'^iP)  ■■■  t\{P-Q] 

r,iP-\}  ■■■  r;(l)  x'^(P+\)  •••  x\{P-Q^\] 

r,(2.V-l)  •  ■  r;(2.V-l-P)  x’^(2.V  -  1)  x’^(2.V  -  1  -  C? ) 

1 2.. 

1)  . 

where  r^in)  =  r^^in)  *  and  T#,(n)  =  r(— n)  •  h\{ri)  and  which  ran  be  solved  as 

before. 
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4.  Maximum  Likelihood  Techniques 

Maximum  likelihood  estimation  of  parameters  is  the  statistical  standard 
against  which  the  performance  of  other  estimators  is  measured.  This  estimator  makes 
use  of  all  the  useful  statistical  information  available  in  a  given  set  of  data  [Ref.  45.  p. 
73].  Difficulties  arise,  however,  in  the  implementation  of  the  maximum  likelihood  es¬ 
timator  (MLE).  True  maximum  likelihood  estimation  requires  exact  knowledge  of  the 
conditional  probability  density  function  (PDF)  of  the  observed  data  conditioned  on 
the  parameters  to  be  estimated.  This  conditional  PDF  must  then  be  simultaneously 
maximized  with  respect  to  all  parameters  being  estimated.  In  practice,  most  efforts 
to  employ  maximum  likelihood  estimation  make  simplifying  assumptions  about  the 
nature  of  the  input  data  to  derive  a  useful  algorithm.  Such  techniques  are  usually 
called  approximate  maximum  likelihood  methods. 

The  approximate  MLE  chosen  for  this  thesis  is  due  to  Akauke  [Ref.  31). 
The  brief  developemcnt  of  this  algorithm  provided  below  follows  that  of  Kay  [Ref.  4. 
Ch.  9.10j.  Additional  backround  on  maximum  likelihood  estimation  can  be  found  in 
[Ref.  1,  8.  32,  45]  and  references  therein. 

Given  a  sequence  of  independent  random  variable  observations,  T(n).  and 
a  corresponding  set  of  parameters  to  be  estimated,  0*,  the  desired  set  of  estimates  for 
the  Ot's  is  the  one  for  which  the  observed  data  set  is  the  most  likely.  In  terms  of  the 
conditional  probability  density  or  likelihood  function, 

p(i(0),x(l),---.x(A'  -  -  - 

the  desired  set  of  estimates  is  the  one  which,  for  a  given  set  of  x(n),  is  maximized 
In  the  case  of  pwle-zero  modeling,  «ui  expression  for  the  observed  data  sequence's 
joint  probability  density  function  conditioned  on  the  model  parameters,  a,  and  b,,  is 
required.  To  obtain  such  an  expression,  it  is  generally  assumed  that  the  observed  data 
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is  Gaussian.  If,  further,  it  is  assumed  that  the  input  to  the  pole  zero  model  is  white 
Gaussian  noise  and  that  the  data  record  length  is  much  longer  than  the  transient 
response  due  to  initial  conditions,  the  conditional  PDF  for  the  observed  data  can  be 
arrived  at  fairly  directly. 

The  joint  PDF  for  the  zero  me2in  white  Gaussian  input  sequence.  u(n),  is 
of  the  form 


p(u(0),u(l] 


,=o  yj2ral 


I 


(2.-36) 


where  is  the  variance.  Now  the  density  function  for  •  •  .x(,V  -  1)  con¬ 

ditioned  on  the  tffc’s,  can  be  found  from  (2.36)  through  the  standard  linear  transfor¬ 
mation. 


p(x(0).x(l),  --,x(.V-  D)  =  p(Uy(0).U/(l)..  -.U/(.V  -  1))1J!. 


(2.37) 


where  u/fn)  is  the  inverse  filter  relationship  of  the  original  pole-zero  difference  equa¬ 
tion.  (2.1).  Specifically 

P  1  <3 


1  1  ^ 
Uf(n)  =  r  ~  r 

^  .=0 


(2.,3S) 


and  J  is  the  Jacobian  of  the  linear  transformation  u;.  To  simplify  the  transformation 
assume  that  the  pole- zero  filter  in  (2.38)  has  been  expressed  as  its  equivalent  all-pole 
filter. 


pAt" 

~  -  i) 


•  sO 


I  he  final  result  of  the  linear  transformation  (2.37)  will  then  be 

p(x(0),  x(  1  )■  •  ■  ■ .  x(A'  -  1 ))  =  n*  -^cxp(  -  ^ > 

nVo 

VVhen  (2.40)  is  expressed  in  terms  of  the  pole-zero  parameters. 


!2,3d) 


(2.40) 


-  ■  ■  ,Op,  •  •  •  ,  1>Q. 
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for  the  purpose  of  maximization  the  relationship  is  highly  nonlinear.  Note  that  the 
maximization  of  (2.40)  requires  the  minimization,  over  all  possible  a,'s  and  6*'s,  of 
the  inverse  filter  error,  u/(n). 

Akaike  employed  the  Newton-Raphson  iterative  method  for  minimizing 
(2.40).  This  method  requires  the  computation  of  Gradient  and  the  Hessian  of  (2.38) 
at  each  iteration  to  generate  the  estimate  updates 


Using  frequency  domain  arguments,  Akaike  was  able  to  provide  expressions  for  the 
above  partial  derivatives  in  terms  of  Fourier  transforms  which  can  in  turn  be  expressed 
in  terms  of  linear  filtering  operations. 

Note  that  there  are  several  key  assumptions  made  above  which  must  be 
valid  for  this  method  of  approximate  NILE  to  apply: 

1.  The  data  are  real,  Gaussian,  and  zero  mean. 

2.  The  data  record  is  large.  This  is  to  avoid  end  effects  of  assuming  all  data  values 
outside  the  data  record  are  zero  when  filtering  the  data. 

3  The  poles  and  zeros  are  not  close  to  the  unit  circle.  This  is  to  avoid  long 
transients  due  to  the  initial  conditions  which  are  ignored.  (They  are  a.ssumed 
to  be  known  and  are  set  equal  to  zero.) 

.-\t  first  glance  these  a.ssumptions  would  seem  to  indicate  that  this  method  is  inappro 
pnate  to  transient  modeling.  However,  inverse  filter  error  retains  its  meaning  when 
ronsidering  transient  waveforms;  the  ideal  inverse  filtering  result  for  a  transient  signal 
IS  a  single  impulse  rather  than  the  minimum  variance  random  .sequence  expected  for 
a  stochastic  process.  In  fact,  this  is  exactly  the  appproach  taken  by  Jackson  [Ref.  26. 
pp.  276-278]  in  extending  Judell’s  m^Lximum  likelihood  method  [Ref.  .33]  to  impulse 
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response  data.  While  this  extension  seems  rather  ad  hoc,  we  will  find  that  such  ap¬ 
proximate  maximum  likelihood  methods  can  be  effective  in  modeling  transient  signals 
cis  impulse  responses.  A  key  limitation  imposed  by  dealing  with  deterministic  data 
is  that  reliance  on  inverse  filter  error  excludes  signals  that  must  be  modeled  by  non¬ 
minimum  phase  systems.  In  contrast,  the  restrictions  of  long  data  record  and  weak 
poles  and  zeros  (not  near  the  unit  circle)  no  longer  apply.  Data  records  end  effects 
and  initial  condition  transient  effects  should  have  no  impact  since  the  assumption  of 
zero  valued  data  outside  the  range  of  data  is  correct  if  the  data  is  chosen  to  end  after 
most  of  the  energy  of  the  transient  is  dissipated. 

C.  IMPLEMENTATION 

.•Ml  modeling  algorithms  were  implemented  using  the  interactive  language  PRO 
.M.ATLAB  from  The  Mathworks.  Inc.  on  SUN  workstations  except  for  the  .-Xkaike 
.MLE  algorithm  which  wa.s  implemented  in  FORTR.^.N  using  a  program  adapted 
from  [Ref.  4.  Ch.  10).  The  FORTR.W  program  was  also  implemented  on  a  SUN 
workstation  with  a  FORTRAN  77  compiler  All  graphics  were  generated  in  PRO¬ 
MAT  LAB. 
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III.  MODELING  PERFORMANCE 

A.  PERFORMANCE  CONSIDERATIONS 

Each  of  the  pole-zero  modeling  techniques  presented  in  Chapter  II  is  effective 
when  modeling  a  signal  that  is  truly  the  impulse  response  of  a  pole- zero  system  with 
no  noise  present  and  with  the  system  order  known  However,  real  world  transient 
data  rarely  possess  such  characteristics.  Real  world  signals  of  all  types  are  notoriously 
uncooperative  in  fitting  the  signal  models  proposed  to  describe  them  Reasons  that 
this  may  be  true  for  transient  data  include: 

1.  Inappropriate  selection  of  model  type  or  modeling  algorithm. 

•  Linear  versus  non-linear  models. 

•  .Minimum  phase  versus  non-minimum  phase  rational  models. 

2.  The  transient  is  time  shifted  because  of  and  inappropriate  selection  of  the  data 
record  starting  point  due  to  the  presence  of  noise. 

■5.  The  assumption  of  impulse  system  excitation  is  a  poor  approximation. 

•1.  Incorrect  selection  of  model  order. 

•').  .Noi.se  is  present  in  the  signal. 

The  test  sequences  used  in  this  chapter  are  all  obtained  as  the  impulse  response  of 
linear  pole-zero  systems,  therefore,  the  question  of  linear  versus  non-linear  model  typie 
will  not  be  at  issue.  For  the  other  problems,  effective  transient  modeling  requires  both 
selecting  the  appropriate  algorithm  and  understanding  how  to  use  that  algorithm  to 
its  greatest  advantage.  The  next  section  describes  the  test  sequences  u.sed  in  this 
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chapter.  Subsequent  sections  address  how  diflRculties  encountered  in  the  pole-zero 
modeling  of  impulse  response  data  can  arise  and  how  they  can  be  overcome. 

The  effects  of  data  selection,  non-minimum  phase  systems,  non-impulse  excita¬ 
tion.  and  incorrect  model  order  on  modeling  will  initially  be  considered  for  signals 
observed  with  no  added  noise  present.  The  effect  of  modeling  a  signal  in  which  addi¬ 
tive  noise  is  present  is  considered  separately.  We  will  see  that  situations  which  modify 
a  linear  pole-zero  system  often  lead  to  another  pole-zero  system.  This  new  system 
usually  has  the  same  number  of  poles  in  the  same  locations  but  with  different  and 
possibly  additional  zeros. 

B.  TEST  SEQUENCES 

The  test  impulse  response  sequences  are  generated  using  pole-zero  models  taken 
from  Kay  [Ref.  4j.  The  test  sequence  ARM.\3  uses  one  of  Kay's  models  directly  while 
the  test  sequences  ARM.A.4  LF.  .\RM.\3  N’.M.  and  .■\RM.A4  CL  are  from  Kay  models 
which  have  betui  modified  to  enhance  (he  illustration  of  certain  points.  The  unit 
impulse  response  and  pole-zero  plot  of  each  test  sequence  model  is  shown  in  Figures 
3.1a-h  The  model  coefficients  for  these  sequences  are  listed  in  Table  3.1. 


Model  i 

Model  Coefficients  (oq  = 

to  —  1  1 

1 

!  <*i 

lOz 

03 

0, 

t, 

ft-: 

AR.MA3  I  -2.760 

-2.654 

0.924 

-0.900 

O.SIO  ; 

ARMA4  LF  i  -3.035 

4.002 

-2.727 

0.778 

-0.200 

0.040 

ARM  A3  N.M  1  -2.760 

,  I 

3.809 

-2.6.54 

0.924 

-4.tsl^ 

25.000  ' 

AR.MA4  f’L  1  -2.67S 

3.700 

-2.56,34 

0.917 

-0.200 

0.040 

TABLE  3.1:  Table  of  test  sequence  coefficients. 
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Figure  3.1:  The  test  sequence  ARMAS,  (a)  Impulse  response  plot  and  (b) 
pole-zero  plot. 
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Figure  3.1:  continued  The  test  sequence  ARMA4  LF.  (c)  Impulse  response 
plot  and  (d)  pole-zero  plot. 
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Figure  3.1:  continued  The  te«t  sequence  ARMAS  NM. 
sponse  plot  and  (f)  pole-zero  plot. 
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Figure  3.1:  continued  The  test  sequence  ARMA4  CL.  (g)  Impulse  response 
plot  and  (h)  pole-zero  plot. 
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C.  ESTIMATING  NUMERATOR  COEFFICIENTS 


1.  Data  Selection — Time  Shifts  and  Initial  Conditions 

Defining  the  range  of  data  to  be  used  in  modeling  is  an  important  and 
usually  staightforward  exercise.  The  assumption  that  a  signal  represents  the  impulse 
response  of  a  linear  pole-zero  system,  however,  implies  some  very  specific  properti<-s 
about  the  initial  few  points  of  that  signal.  Each  method  of  modeling  the  transfer 
function  numerator  coefficients  in  Chapter  II  reacts  differently  when  the  beginning 
points  of  the  impulse  response  being  modeled  are  degraded.  Because  real  worhi 
transients  do  not  usually  exhibit  the  instantaneous  ri.se  time  of  an  ideal  imjuilsf* 
response  and  because  noise  is  usually  present,  choosing  the  precise  data  range  for  a 
transient  such  as  that  illustrated  in  F'igure  3.2  is  often  a  very  uncertain  task.  I  wo 
possible  outcomes  when  the  transient  starting  point  is  chosen  incorrectly  are: 

1.  I  he  starting  point  is  chosen  before  the  signal  begins  so  that  earl>’  data  %alues 
are  unri'lated  land  presumably  of  lower  amplitude i  to  the  imjiulse  response  to 
be  matched  (e.g.  these  points  may  consist  of  noise i, 

'I  lie  starting  point  is  chosen  late  in  which  ca.s<-  early  valut's  of  tlie  impuKe 
re-^ponse  are  lost. 

In  the  first  ra.se.  an  adequate  number  of  additional  model  zero^  <  an  ai  count 
for  the  delay  in  the  impulse  response  .Assuming  that  the  startinc  point  for  the  data  is 
rea.sonably  close  to  the  true  beginning  of  the  impulse  response,  an\'  spectral  features 
introduced  by  the  unrelated  early  data  points  will  probably  not  have  enough  energv  to 
significantly  alter  the  spectrum  of  the  impulse  response.  I’nder  these  circumstances 
all  methods  can  effectively  find  the  jioles.  It  is  important  to  note,  however,  that  if 
an  insufficient  numfier  of  zeros  to  account  for  the  imposed  rlelay  is  used,  then  some 
of  the  equations  that  are  generated  in  Prony's  method  become  invalid.  W  hen  these 


.30 


-2000^ - - - 

0  500  1000 

n 

Figure  3.2;  An  example  of  a  laboratory  generated  acoustic  transient.  Note 
the  difficulty  in  determining  a  precise  starting  point  for  the  transient. 


equations  are  solved,  the  invalid  equations  can  drastically  degrade  pole  estimates.  To 
see  this  we  cm  apply  Prony’s  method  to  a  system  with  the  true  orders  P  =  4  and 
Q  =  2.  Assume  the  signal  is  delayed  by  inserting  three  zeros  at  the  oeginning  of  the 
data  so  that  the  original  point  i{0)  is  now  the  fourth  data  point.  If  we  choose  P  =  4 
and  Q  =  3  in  constructing  the  data  matrix  of  (2.2),  the  resulting  set  of  equations  will 
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Note  that  when  the  lower  partition  is  taken  to  find  the  ai,  coefficients,  the  first  two 
equations  of  the  lower  partition  are  invalid.  Compounding  the  problem  is  that  the 
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invalid  equations  occur  in  the  high  energy  portion  of  a  transient  signal.  To  overcome 
this  effect  a  numerator  order  of  at  le^lst  five  is  required. 


In  finding  the  numerator  coefficients  of  a  delayed  impulse  response  one 
of  three  outcomes  is  possible  depending  on  the  estimation  technique  chosen,  first, 
any  technique  which  depends  directly  on  the  initial  values  of  the  data  sequence  (for 
example  (2.5))  will  be  ineffective.  Second,  methods  which  rely  on  the  autocorrelation 
of  the  residual  sequence  will  result  in  approxirnatelv  the  true  zeros  of  the  system  under 
study.  The  original  time  series  will  not  be  matched  directly.  This  case  is  illustrated 
in  Figure  3. .3.  The  resulting  impulse  response  is  an  undelayed  version  of  the  signal. 
Finally,  signal  matching  techniques,  iterative  prefiltering  and  Shank’s  method,  result 
in  zeros  not  related  to  the  original  undelayed  model  but  which  provide  the  best  o\’eral! 
time  domain  match  of  the  delayed  signal.  This  is  shown  in  Figure  3.4. 

The  ca,se  of  choosing  the  data  record  to  far  to  the  right  and  thus  trun¬ 
cating  the  first  points  of  an  impulse  response  will  again  have  little  effect  on  pole 
estimation-  This  situation  corresponds  the  same  system  with  (non-zero)  initial  con¬ 
ditions  imposed.  Since  initial  conditions  are  acounted  for  in  the  numerator,  the  zeros 
are  significantly  altered.  Here  the  previous  discussion  regarding  finding  the  under¬ 
lying  mode!  zeros  versus  obtaining  a  good  match  in  the  time  domain  still  applies 
with  one  exception:  the  direct  calculation  of  the  6*  coefficients  from  i2.5i  will  now 
he  effective. 

2.  Non-minimum  Phase  Modeling 

In  [Hef.  4r)]  it  is  demonstrated  that  the  appropriate  discrete  model  tif 
a  sampled  analog  waveform  is  often  represented  by  a  transfer  function  with  zeros 
outside  the  unit  circle.  Many  of  the  techniques  that  are  currently  available  are  pur 
[)f>sefully  structured  to  rltmtnair  uch  non-minimum  phase  models.  In  power  spectrum 
estimation,  a  minimum  pha.se  system  with  the  same  frequency  response  magnitude 
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as  a  non-minimum  pha^e  system  results  in  the  same  estimated  spectrum.  Thus  for 
spectrum  estimation  an  equivalent  minimum  phase  system  is  satisfactory.  In  fact, 
the  assumption  underlying  stoch«tic  modeling  techniques,  namely  white  Gaussian 
noise  input,  guarantees  that  all  transfer  function  combinations  of  minimum  phase 
and  maximum  phase  zeros  are  equivalent.  By  convention,  stochastic  modeling  tech¬ 
niques  always  choose  the  minimum  phase  model  so  that  the  important  statistic  of 
inverse  filler  error  is  available.  In  time  domain  based  applications,  however,  incor¬ 
rectly  choosing  the  model  phase  caji  seriously  degrade  system  performance.  Fields 
such  as  seismic  deconvolution,  channel  equilization.  control,  and  matched  filter  de¬ 
sign.  generaily  require  identification  of  the  correct  model  phase  [Ref.  47.  48.  491. 
.Mso.  we  will  see  in  Chapter  I\'  that  effective  modeling  of  real  world  acoustic  signals 
frequently  requires  non-minimum  phase  models. 

A  number  of  modifications  to  the  basic  stochastic  model  have  been  in¬ 
troduced  to  allow  selection  of  the  model  with  the  correct  phase.  These  techniques 
generally  involve  changing  the  Gaussian  nature  of  the  input  noise  [Ref.  .dU]  and  often 
employ  higher  than  second  order  moments  or  cumulants  [Ref.  29.  30j.  However,  when 
a  model's  impulse  response  has  effectively  matched  a  signal  in  the  time  domain,  the 
resulting  phase  is  immediately  known  to  be  the  correct.  .Allowing  for  the  possibilitx 
of  non-minimum  pha.se  models  places  significant  limitations  on  the  modeling  methods 
which  may  be  used.  Techniques  which  rely  on  inverse  filtering  I  maximum  liklihood 
methods )  are  not  applicable  since  the  inverse  of  a  non-minimum  phase  system  is  unsta¬ 
ble.  .Also,  techniques  which  use  correlation  information  to  calculate  the  coefficients 
(spectral  factorization  and  Durbin's  method)  will  give  poor  results  since  correlation 
data  does  not  preserve  pha.se  information.  Figure  .3..')  shows  the  difficulty  encountered 
when  Durbin’s  method  attempts  to  model  a  non-minimum  phase  system.  The  best 
Durbin’s  method  can  do  is  produce  the  spectrally  equivalent  minimum  phase  version 
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of  a  maximum  phtise  system  since  the  autoregressive  modeling  techniques  on  which 
it  relys  can  only  produce  zeros  within  the  unit  circle.  In  contrast.  Figure  3.6  shows 
that  Shank’s  method  is  able  to  find  the  correct  model. 

3.  Numerator  Modeling  Summary 

Table  3.2  provides  a  brief  summary  of  the  modeling  properties  of  the  nu¬ 
merator  coefficient  modeling  techniques  considered  in  this  thesis.  Since  the  goal  in 
Chapter  IV  is  to  perform  time  domain  modeling  of  the  acoustic  transients  being 
considered,  those  methods  which  provide  the  best  time  domain  match  between  the 
original  signal  and  the  model  impulse  response  are  preferred. 
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TABLE  3.2:  Summary  of  the  capabilities  and  limitations  of  numerator 
modeling  methods. 


D.  NON-IMPULSE  EXCITATION 

In  any  real  world  system  the  assumption  of  a  unit  impulse  input  is  approxi 
mate.  If  the  duration  of  the  excitation  waveform  is  small  relative  to  the  period  of  the 
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Figure  3.3;  The  right  shifted  sequence  ARMA3,  poles  modeled  using 
LSMYWE  and  zeros  modeled  using  Durbin’s  method,  (a)  Time  signal 
plot  and  (b)  pole-zero  plot.  Durbin’s  method  does  not  account  for  the 
time  delay  but  instead  finds  the  underlying  system’s  true  zeros. 
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Figure  3.4:  The  right  shifted  sequence  ARMAS,  poles  modeled  using 
LSMYWE  and  zeros  modeled  using  Shank’s  method,  (a)  Time  signal 
plot  and  (b)  pole-zero  plot.  The  least  squares  method  does  not  find  the 
underlying  system’s  true  zeros  but  rather  achieves  the  best  overall  time 
domain  match. 
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Figure  3.5:  The  non-minimum  phase  sequence  ARMAS  NM,  poles  mod¬ 
eled  using  LSMYWE  and  zeros  modeled  using  Durbin's  method,  (a)  Time 
signal  plot  and  (b)  pole-zero  plot.  Durbin’s  method  cannot  model  zeros 
outside  the  unit  circle. 
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Figure  3.6.  The  non-minimum  phase  sequence  ARMAS  NM,  poles  mod¬ 
eled  using  LSMYWE  and  xeros  modeled  using  Shank’s  method,  (a)  Time 
signal  plot  and  (b)  pole-zero  plot.  The  least  squares  method  is  effective 
at  modeling  zeros  outside  the  unit  circle. 
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lowest  oscillation  frequencies  present  then  the  assumption  is  justified.  However,  it  is 
reasonable  to  expect  this  criterion  will  often  not  be  satisfied.  Non-impulse  excitations 
which  may  be  encountered  include: 


1 .  A  long  duration  baseband  type  pulse. 

2.  An  uncorrelated  random  train  of  impulses.  This  model  is  often  used  to  account 
for  reverberation  (echoes)  in  siesmic  deconvolution  [Ref.  51]. 

.3.  A  frequency  swept  input  that  sweeps  through  the  natural  frequencies  of  a  sys¬ 
tem.  This  model  is  usually  considered  in  conjunction  with  the  starting  and 
stopping  of  rotating  machinery. 

If  the  system  input  were  known,  the  modeling  problem  could  be  formulated  as  a  sys¬ 
tem  identification  problem.  When  no  information  about  the  system  input  is  available, 
other  means  must  be  found  to  deal  with  this  problem. 

1.  Baseband  Pulse  Excitation 

The  effect  of  modeling  a  transient  signal  from  a  linear  system  in  which  the 
input  is  a  long  duration,  baseband-type  pulse  can  best  be  understood  as  filtering  by  a 
finite  impulse  response  (FIR,  all-zero)  filter  as  illustrated  in  Figure  3.7  The  spectral 
properties  of  the  original  time  series  are  windowed  by  the  frequency  response  of  the 
FIR  filter  coefficients.  For  this  type  of  pulse,  the  effect  is  that  of  low  pass  filtering. 
Thus  high  frequencies  are  attenuated  relative  to  low  frequencies.  If  no  frequency 
component  exists  below  the  FIR  filter's  cutoff  frequency  then  the  original  spectrum 
IS  altered  according  to  the  side  lobe  structure  of  the  FIR  filter. 

The  new  model  that  results  can  be  viewed  a.s  a  system  with  the  original 
model  poles  but  with  new  numerator  polynomial  coefficients  that  are  the  result  cf 
convolving  the  FIR  filter  coefficients  with  the  original  numerator  polynomial  coeffi 
cients.  This  results  in  a  higher  order  polynomial,  hence  more  zeros  than  were  present 
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Figure  3.7:  One  way  to  view  a  linear  system  excited  by  a  baseband  pulse. 


in  the  original  model  will  be  required.  Observe  how  the  original  impulse  response  of 
Figure  3.7a  is  altered  to  Figure  3.7b  when  a  nine  element  triangular  excitation  pulse 
is  used.  The  dotted  lines  in  Figures  3.7b.c.d  indicate  the  modeling  results  obtained 
when  none,  two.  and  four  extra  zeros,  respectively,  are  used  in  the  estimated  pole-zero 
model.  In  this  case  four  extra  zeros  prove  sufficient  to  account  for  the  input  pulse. 
The  corresponding  model  spectrum.  Figure  3.7e  illustrates  the  attenuation  caused  by 
the  bziseband  excitation.  The  pole- zero  plot  in  Figure  3.7f  shows  that  non-minimum 
phase  zeros  were  required  to  achieve  an  effective  time  domain  match. 

2.  Random  Impulse  Train  Excitation 

If  the  input  to  a  linear  system  is  an  uncorrelated  random  tram  of  impulses 
then,  although  the  time  series  may  be  significantly  different  from  the  original  impulse 
response,  the  autocorrelation  function  of  the  signal  is  theoretirall\’  unaltered  except 
for  a  scaling  factor.  This  is  because  the  autocorrelation  of  an  uncorrelated  impulse 
train  is  a  scaled  unit  impulse.  Thus  modeling  methods  which  rely  on  correlation 
information  should  be  effective.  However,  over  a  finite  time  interval  it  is  unlikely  that 
a  random  sequence  will  be  truly  uncorrelated.  As  the  impulse  train  becomes  correlated 
the  situation  will  be  equivalent  to  the  baseband  pulse  ca.se  described  afiove. 

3.  Frequency  Swept  Excitation 

The  output  of  a  system  excited  with  a  frequency  swept  signal  depends  on 
the  rate  at  which  the  sweep  occurs.  A  slow  sweep  will  result  in  a  series  of  transient 
events,  each  at  a  specific  resonant  frequency  of  the  system.  These  events  can  each 
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be  modeled  separately.  When  the  sweep  rate  is  rapid,  all  natural  modes  will  appear 
much  as  if  the  input  were  an  impulse  except  that  the  phase  relationships  of  the  various 
components  may  be  changed. 

E.  MODEL  ORDER  SELECTION 

In  studies  of  rational  modeling  the  issue  which  continues  to  be  the  most  con¬ 
founding  is  that  of  model  order  selection.  The  proposed  methods  which  have  a  sound 
theoretical  basis  (e.g.  [Ref.  52,  53.  54])  are  very  difficult  to  actually  implement. 
These  methods  are  invariably  related  to  maximum  likelihood  concepts  and  therefore 
rely  heavily  on  inverse  filter  statistics.  This  implies  that  for  model  order  evaluation 
the  inverse  filter  error  must  be  calculated  over  all  possible  model  orders.  Then  the 
model  order  and  inverse  filter  error  which  minimize  .some  function  of  the  two  is  se¬ 
lected.  The  case  of  non-minimum  phase  systems  is  even  more  intractable  since  the 
inserse  filter  of  such  systems  is  unstable. 

A  more  attractive  but  less  understood  method  is  to  initiall\’  overdetermine  a 
system  and  then  allow  the  modeling  algorithm  to  indicate  the  correct  model  order. 
.A  method  proposed  by  (’adzow  [Ref.  39]  uses  singular  value  decomposition  to  aid  in 
model  order  selection  in  the  denominator  of  all-pole  and  pole- zero  models.  Kumeresan 
and  Tufts  [Ref.  21]  have  shown  that  when  Prony’s  method  is  applied  to  exponentially 
damped  sinusoids  reversed  in  time,  valid  poles  occur  outside  the  unit  circle  and  excess 
poles  occur  inside  the  unit  circle.  In  both  of  these  methods  the  denominator  order 
IS  initially  overdetermiried  and  the  modeling  method  then  provides  the  corre»ct  order. 
For  this  thesis  both  of  these  methods  were  applied  to  a  number  of  transients.  They 
were  effective  when  applied  to  the  noiseless  impulse  responses  of  true  pole- zero  svstems 
but  were  not  robust  in  the  presence  of  noise  or  when  many  narrowband  components 
are  present.  There  is  no  similar  guidance  for  determining  the  numerator  model  order. 
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Figure  3.8:  The  sequence  ARMA4  LF  with  (a)  unit  impulse  excitation 
and  (b)  triangular  pulse  input  modeled  with  correct  model  orders  P  =  4 
and  Q  =  2  using  Prony’s  method  and  Shank’s  method. 
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Figure  3.8:  continued  The  sequence  ARMA4  LF  with  triangular  pulse 
input  (c)  modeled  with  Prony’s  method  and  least  squares  identification 
with  orders  P  =  4  and  Q  =  4  (d)  modeled  with  Prony’s  method  and  Shank’s 
method  orders  P  =  4  imd  Q  =  6. 
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Figure  3.8:  continued  The  sequence  ARMA4  LF  with  triangular  input  (e) 
the  original  impulse  response  and  triangular  output  model,  order  (4,6), 
spectrums  and  (f)  the  corresponding  poIe-*ero  plot.  Note  how  the  base¬ 
band  pulse  has  attenuated  the  higher  frequencies. 


Perhaps  most  desirable  would  be  a  modeling  technique  that,  when  overdeter¬ 
mined,  caused  excess  poles  and  zeros  to  either  cancel  or  migrate  well  away  from  the 
unit  circle  (all  the  way  to  the  origin  ideally).  Tummala  [Ref.  23]  has  demonstrated 
some  success  with  this  concept  using  an  iterative  algorithm  to  solve  the  least  squares 
identification  problem.  Observing  the  degree  to  which  the  various  transient  modeling 
techniques  of  this  thesis  exhibit  this  behavior  is  the  approach  taken  in  the  following 
comparison.  This  behavior  was  observed  by  modeling  the  ARMA3  test  sequence  with 
different  combinations  of  numerator  and  denominator  order.  Table  3.3  summarizes 
the  results  obtained.  A  ‘Y’  in  Table  3.3  indicates  a  modeling  technique  achived  an  ex¬ 
act  time  domain  match  between  the  model  impulse  respon.se  the  modeled  signal.  The 
most  notable  negative  result  is  that  both  Prony's  method  and  LSMVVVE  were  ineffec¬ 
tive  when  the  numerator  and  denominator  were  overdetermined.  .Mso.  excess  poles 
without  enough  zeros  causes  problems  for  correlation  domain  iterative  prefiltering. 
Figure  3.9  is  an  example  of  the  results  obtained  using  LSMYWE  when  excess  poles 
and  zeros  are  present.  However,  when  these  results  were  used  to  initialize  iterative 
prefiltering  the  excess  poles  and  zeros  were  handled  effectively. 

The  above  results  and  experience  gained  during  extensive  modeling  trials  lead  to 
the  following  recommended  strateg%'  when  modeling  complex  transients  about  which 
very  little  is  known: 

I.  Be  conservative  in  estimating  the  denominator  order.  Signals  composed  of  nar 
rowband  components  will  generally  be  dominated  by  the  few  components  high¬ 
est  in  energ}'.  Even  if  excess  zeros  are  required  to  deal  with  a  noisy  signal  (see 
the  next  section)  a  small  number  of  excess  zeros  will  usually  suffice  when  us¬ 
ing  LSMYWE.  Pole  detection  techniques  such  as  those  of  Cadzow[Ref.  5],  and 
Kumeresan.  and  Tufts  [Ref.  21]  may  be  helpful. 
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TABLE  3.3;  The  effectiveness  of  thesis  modeling  methods  on  ARMA3 
with  varying  overdetermination  of  model  order.  ‘Y’  indicates  an  exact 
time  domain  match  was  achieved  and  ‘N’  indicates  a  poor  time  domain 
match  resulted. 

2.  Be  expansive  with  estimates  of  numerator  order.  (Iterative  prefiltering  may  be 
an  exception  to  this.  See  Chapter  I\  .)  .Additional  zeros  can  help  compensate 
for  mistakes  made  in  data  selection  and  effects  of  non-impulse  excitation.  I  ri- 
needed  zeros  will  usually  migrate  away  from  the  unit  circle.  Note  that  when 
using  Prony  s  method  or  LS.MYWE.  one  numerator  order  can  be  used  when 
calculating  the  denominator  coefficients  and  amther  numerator  order  can  be 
chosen  when  finding  the  numerator  coefficients.  Using  Shank's  method  with 
Prony ’s  method  or  LSMYWE  provides  the  most  flexibility  when  a  model  is 
attempting  to  account  for  erroneous  assumptions  (e.g.  data  selection)  which 
may  have  been  made  by  the  user.  The  above  recomendations  are  primarilv 
intended  for  Prony ’s  method  and  LSMYWE.  However,  using  any  iterative  tech¬ 
nique  requires  reasonable  initial  estimates  which  are  usually  arrived  at  using 
linear  estimation  techniques. 
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Figure  3.9:  PoIe**ero  plots  for  the  sequence  ARMAS  modeled  with  overde¬ 
termined  model  orders  of  P=10  and  Qs=10.  (a)  Ismywe  and  Durbin's 

method  results  in  an  unstable  system,  (b)  iterative  preflltering  provides 
an  effective  estimation  of  the  correct  poles  and  zeros  and  an  excellent  time 
domain  match  (not  shown). 
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F.  NOISE  PERFORMANCE 


Kay  [Ref.  55]  has  shown  that  for  an  all-pole  process,  additive  white  noise  h£is 
the  effect  of  introducing  zeros  which  migrate  from  the  origin  to  the  model  poles  as 
the  signal-to-noise  ratio  is  decreased.  Alternatively,  if  zeros  are  not  introduced  into 
the  model,  poles  move  toward  the  origin  as  the  signal-to-noise  ratio  is  reduced  [Ref. 
56].  In  either  case  the  overall  effect  is  a  loss  of  spectral  resolution.  This  result  extends 
directly  to  pole-zero  processes  in  white  noise.  Therefore  one  important  feature  of  any 
modeling  technique  is  it’s  ability  to  resolve  spectral  components  in  the  presence  of 
noise. 

1.  Rational  Modeling  of  Noisy  Data 

Several  authors  have  recently  addressed  the  difficulties  associated  with 
modeling  a  time  series  in  which  additive  white  noise  is  present.  Kay  [Ref.  57]  has 
noted  that  when  the  variance  of  the  white  noise  can  be  estimated,  its  effect  can 
be  removed  from  the  zeroth  autocorrelation  lag  to  improve  the  resolution  of  pole 
frequencies  in  correlation  based  autoregressive  modeling.  .Alternatively,  the  zeroth 
autocorrelation  lag  may  be  avoided  by  using  (2.18).  sometimes  called  the  htgh  ordrr 
Yule- Walker  Equations,  for  Q  >  P  or  by  eliminating  the  rows  containing  the  zeroth 
lag  in  certain  lea.st  squares  formulations  [Ref.  44].  If  the  resulting  high  order  A'ule- 
Walker  equations  yield  a  matrix  that  is  well  conditioned,  this  procedure  is  equivalent 
to  that  of  Kay  above  [Ref.  58).  Even  with  good  correlation  estimates,  however,  the 
matrix  Rx  is  not  guaranteed  to  be  positive  definite  (i.e.,  invertible).  Consequently, 
Ra  is  more  likely  to  be  a  poorly  conditioned  matrix.  This  is  not  a  problem  in  the 
least  squares  formulation.  (2.19),  since  the  very  act  of  choosing  such  a  formulation 
implies  we  expect  a  solution  that  is  approximate. 

A  number  of  authors  have  noted  that  overdetermining  the  number  of  model 
poles  can  have  a  profound  effect  when  modeling  noisy  signals  [Ref.  55,  56,  21, 36,  39]. 
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The  explanation  for  this  is  that  the  extra  poles  eire  able  to  model  the  flat  noise 
spectrum  of  the  signad,  preventing  that  portion  of  the  signal  from  interfering  with 
the  poles  needed  for  narrowband  modeling.  Finally,  singular  value  decomposition  has 
been  used  to  aid  resolution  by  overcoming  the  conditioning  problems  that  can  result 
from  pole  overdetermination  and  high  order  Yule- Walker  equations  [Ref.  5.  21.  36, 
59,  60]. 

Little  work  has  appeared  which  analyzes  or  demonstrates  the  effectiveness 
of  iterative  prefiltering  in  the  presence  of  additive  white  noise.  Stoica  and  Soder- 
strom  [Ref.  61]  have  shown  that  iterative  prefiltering  will  perform  well  for  very  long 
stochiistic  data  sequences  in  the  presence  of  additive  white  noise  if  reasonable  initial 
estimates  are  available  to  start  the  iterations.  They  also  show  that  for  arbitrary  initial 
estimates  iterative  prefiltering  is  not  guairanteed  to  converge.  Tufts  and  Kumeresan 
[Ref.  62]  have  demonstrated  superior  performance  for  approximate  MLE  over  linear 
prediction  methods  when  modeling  sinusoids  in  white  noise.  Using  the  approximate 
MLE  method  of  Box  and  Jenkins.  Cadzow  found  that  such  a  method  performed  com 
parably  to  the  overdetermined  high  order  V'ule- Walker  equations  for  the  sum  of  two 
second  order  all-pole  processes  in  white  noise  but  with  higher  variance. 

The  above  work  suggests  numerous  possibilities  in  dealing  with  noisy  im¬ 
pulse  response  data.  The  effectiveness  of  these  techniques  is  evaluated  empirically  in 
the  following  section  by  modeling  noisy  test  sequences. 

2.  Discussion  and  Examples 

The  examples  used  to  illustrate  modeling  performance  are  summarized 
in  Table  3.4.  All  noise  modeling  examples  were  conducted  using  the  test  sequence 
ARMA4  CL  (see  Table  3.1).  Figure  3.10  shows  that  Prony’s  method  has  difficulty 
even  when  the  SNR  is  as  high  as  30  dB  although  the  excess  poles  introduced  for 
Figure  3.10b  dramatically  increase  modeling  effectiveness.  The  moderate  impact  of 
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insuring  that  Q  >  P  is  illustrated  for  LSMYWE  in  Figure  3.11  in  that  the  higher 
frequency  pole  moves  out  toward  the  unit  circle  when  the  modeling  order  is  altered 
from  (4,2)  to  (4,4).  As  with  Prony’s  method,  overestimating  the  number  of  poles 
dramatically  improves  modeling  in  Figures  3.12  and  3.13.  In  Figure  3.14.  the  two 
previous  noise  modeling  examples  are  again  modeled  with  the  smallest  singular  value 
set  to  zero  to  yield  a  data  matrix  of  rank  P.  The  technique  is  highly  effective  in 
both  cases.  The  effectiveness  of  singular  value  decomposition  when  overdetermining 
the  number  of  poles  is  illustrated  by  the  examples  in  Figure  3.15.  This  is  identical 
to  the  result  described  in  [Ref.  21].  Discarding  the  excess  singular  values  aids  both 
in  resolving  the  narrowband  components  present  and  in  reducing  the  variance  of  the 
excess  poles.  Figures  3.16  and  3.17  show  the  haz^u•ds  that  are  possible  when  using 
singular  value  decomposition.  .Any  number  of  singular  values  discarded  other  than 
the  those  necessary  to  achieve  a  data  matrix  rank  equal  to  the  transfer  function 
denominator  order  has  undesirable  side  effects.  Discarding  too  few  singular  values 
such  as  in  the  Figure  3.16  examples  will  resolve  the  desired  narrowband  components 
but  also  lead  to  excess  poles  near  or  even  beyond  the  unit  circle.  This  result  will 
at  best  increase  the  likelihood  of  spurious  narrowband  components  and  can.  in  the 
worst  case,  lead  to  an  unstable  model.  If  too  many  singular  values  are  discarded 
as  in  Figure  3.17b,  the  accuracy  of  narrowband  component  estimates  will  be  badlv 
degraded.  Therefore,  singular  value  decomposition,  w’hile  potentially  very  useful  in 
modeling  transient  signals  in  noise,  must  be  used  with  caution  when  the  true  model 
order  is  unknown  (which  is  basically  always). 

The  iterative  prefiltering,  correlation  domain  iterative  prefiltering,  and 
Akaike  MLE  algorithms  were  initialized  with  several  of  the  preceeding  examples  and 
the  results  are  shown  in  and  in  Figures  3.18  and  3.19.  Iterative  prefiltering  dramati¬ 
cally  improved  the  narrowband  modeling  of  each  example  with  which  it  was  initialized. 
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Figure  3.18a  indicates  that  excess  poles  may  cause  problems  by  being  close  to  the  unit 
circle.  However,  the  ability  of  iterative  prefiltering  to  cancel  excess  zeros  with  excess 
poles  demonstrated  in  this  chapter  may  mitigate  this  effect.  Correlation  domain  iter¬ 
ative  prefiltering  also  improved  the  resolution  of  each  noise  example  on  which  it  was 
tried  although  its  convergence  was  slower  than  and  the  bias  of  its  results  were  higher 
than  standard  iterative  prefiltering.  Also,  correlation  domain  iterative  prefiltering 
showed  a  alarming  tendency  to  place  excess  poles  near  the  true  poles  on  noisy  sig¬ 
nals.  We  experienced  a  great  deal  of  difficulty  in  trying  to  model  noisy  signals  with 
the  Akaike  MLE  algorithm.  The  primary  problem  encountered  was  early  algorithm 
termination  due  to  an  iteration  which  produced  a  non-minimum  phase  zero.  Akaike 
.MLE  was  able  to  improve  the  quality  of  a  ‘good'  Prony  estimate  in  Figure  3.19b. 
Finally.  Figure  3.20  shows  a  pole-zero  modeling  result  corresponding  to  the  example 
in  Figure  3.17a.  Figure  3.20a  shows  the  noisy  signal  used  to  during  modeling  and 
Figure  3.20b  compares  the  resulting  model  to  the  original  noiseless  impulse  response. 

G.  MODELING  PERFORMANCE  SUMMARY 

In  modeling,  the  more  that  is  known  about  the  signal  to  be  modeled,  the  better 
the  choice  of  model  structure  and  modeling  algorithm  that  can  be  made,  The  focus 
of  this  thesis  is  the  modeling  of  real  world  transient  signals  about  which  very  little  is 
known  and  whose  characteristics  can  vary  widely.  The  type  of  model  is  also  set  by  the 
basic  goal  of  this  thesis,  that  is,  a  rational  linear  model,  I’nder  such  circumstances, 
a  sensible  criterion  in  choosing  a  modeling  algorithm  is  robustness.  For  numerator 
modeling.  Shank's  method  is  particularly  robust  in  the  sense  that  it  provide  the 
numerator  coefficients  which  give  the  minimum  norm,  least  squares  fit  based  on  the 
previously  determined  denominator  coefficients.  Shank's  method  is  insensitive  to 
time  shifted  data  records  and  uncertainties  in  model  order.  Durbin's  method  is  also 
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effective  but  suffers  from  the  limitations  that  only  minimum  phase  models  are  possible 
and  the  resulting  models  cannot  account  for  time  shifts.  Since  the  Akaike  MLE 
algorithm  requires  a  minimum  phase  initial  estimate,  Durbin’s  method  can  be  used 
provide  an  such  an  estimate.  Equation  (2.5)  and  spectral  factorization  are  extremely 
sensitive  to  degraded  signals  (e.g.  noisy  signals)  and  so  will  not  be  considered  further. 

Denominator  modeling  with  Prony’s  method  is  extremely  sensitive  to  even  small 
amounts  of  noise.  This  is  only  slightly  less  true  of  Akaike  MLE.  Also.  Akaike  MLE 
cannot  continue  to  iterate  if  a  non-minimum  phase  model  estimate  is  encountered.  In 
contrast,  LSMYWE  and  iterative  prefiltering  are  quite  robust  in  noise  and  iterative 
prefiltering  is  particularly  robust  to  model  order  overdetermination.  Our  conclusion 
is  that  for  modeling  the  real  world  acoustic  transients  of  the  next  chapter.  LSM'^’WE 
with  Shank's  method  and  iterative  prefiltering  are  likely  to  perform  the  best  on  the 
acoustic  transients  considered  in  the  next  chapter. 


Method(order  P,Q) 
(initialization) 
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in  dB 

Bias  Variance  of  Pole 

(true  minus  model)  Estimates 

Figure 

True  Poles 

.6247  ±  .7509> 
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30 
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TABLE  3.4:  Examples  of  pole-zero  modeling  performance  of  the  impulse 
response  test  sequence  ARMA4  CL  in  added  noise. 


X  true  poles 
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Figure  3.10:  Model  pole  scatter  plots  of  twenty  trials  illustrating  the  ef¬ 
fectiveness  of  overdetermining  the  number  of  model  poles  using  Prony's 
method.  Test  sequence  is  ARMA4  CL  with  SO  dB  of  noise,  (a)  Using  the 
correct  model  order,  (4,2)  and  (b)  Using  two  excess  poles,  model  order 
(6,2). 
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Figure  3.11:  Model  pole  scatter  plots  of  twenty  trials  illustrating  the  mod¬ 
erate  benefits  of  zeroth  lag  cc^rrelation  compensation  by  choosing  Q  =  P 
when  using  LSMYWE.  Test  sequence  is  ARMA4  CL  with  SO  dB  of  noise, 
(a)  Using  the  correct  model  order,  (4,2)  and  (b)  model  order  (4,4). 
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X  true  poles 
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Figure  3.12;  Model  pole  icatter  plots  of  twenty  trials  illustrating  the 
dramatic  benefits  of  using  excess  poles  with  LSMYWE.  Test  sequence 
is  ARMA4C.  (a)  Model  order  (6,2)  with  15  dB  of  added  noise  and  (b) 
model  order  (8,8)  with  lOdB  of  added  noise. 
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X  true  poles 
.  estimated  poles 


X  true  poles 
.  estimated  poles 


Figure  3.13;  Model  pole  scatter  plots  of  twenty  trials  illustrating  the 
dramatic  benefits  of  using  excess  poles  with  LSMYWE.  Test  sequence 
is  ARMA4  CL.  (a)  Model  order  (8,8)  with  5  dB  of  added  noise  and  (b) 
model  order  (12,12)  with  5  dB  of  added  noise. 
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Figure  3.14:  Model  pole  scatter  plots  of  twenty  trials  illustrating  the  dra¬ 
matic  benefits  of  setting  the  smallest  singular  value  of  the  data  matrix  to 
zero  (i.e.  adjust  data  matrix  rank  to  P)  for  Prony’s  method  and  LSM^'A\'E. 
Test  sequence  is  ARMA4  CL.  (a)  Prony’s  method  for  model  order  (4,2) 
with  30  dB  of  added  noise  and  (b)  LSMYWE  for  model  order  (4,4)  with 
15  dB  of  added  noise. 
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Figure  3.15:  Model  pole  scatter  plots  of  twenty  trials  illustrating  the  effect 
of  adjusting  the  data  matrix  rank  when  using  excess  poles  for  Prony’s 
method.  Test  sequence  is  ARMA4c.  (a)  Model  order  (8,2)  with  20  dB 
of  added  noise  and  no  rank  acijustment  and  (b)  model  order  (8,2)  with 
20  dB  of  added  noise  with  rank  adjusted  to  =  4  using  singular  value 
decomposition. 
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Figure  3.16:  Model  pole  scatter  plots  of  twenty  trials  illustrating  the  effect 
of  adjusting  the  data  matrix  rank  when  using  excess  poles  for  LSMYWE. 
Test  sequence  is  ARMA4  CL.  (a)  Model  order  (8,8)  with  5  dB  of  added 
noise  with  rank  adjusted  to  +  4  =  8  and  (b)  model  order  (8,8)  with  5 
dB  of  added  noise  with  rank  adjusted  to  +  2  =  6. 
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Figure  3.17:  Model  pole  scatter  plots  of  twenty  trials  illustrating  the  effect 
of  adjusting  the  data  matrix  rank  when  using  excess  poles  for  LSMYWE. 
Test  sequence  is  ARMA4  CL.  (a)  Model  order  (8,8)  with  5  dB  of  added 
noise  with  rank  adjusted  to  P(r»«  =  4  and  (b)  model  order  (8,8)  with  5  dB 
of  added  noise  with  rank  adjusted  to  P,rM  —  1=3. 
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Figure  3.18:  Model  pole  scatter  plots  of  twenty  trials  illustrating  the  abil¬ 
ity  of  iterative  prefiltering  to  improve  the  resolution  of  an  LSMYWE  es¬ 
timate.  In  both  cases  the  iterative  prefiltering  algorithm  was  initialized 
using  LSMYWE  and  Durbin’s  method.  Test  sequence  is  ARMA4  CL.  (a) 
Model  order  (4,4)  with  15  dB  of  added  noise  and  (b)  model  order  (8,8) 
with  5  dB  of  added  noise. 
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X  true  poles 
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Figure  3.19:  Model  pole  ecatter  plot*  of  twenty  trials  on  test  sequence 
ARMA4  CL  illustrating  (a)  the  ability  of  correlation  domain  prefiltering 
to  improve  the  resolution  of  an  LSMYWE  initial  estimate,  order  (8,8),  5 
dB  added  noise  and  (b)  the  ability  of  Akaike  MLE  to  improve  the  resolution 
of  a  Prony’s  method  initial  estimate,  model  order  (6,2),  SO  dB  added  noise. 
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Figure  3.20:  Pole-zero  model  of  the  ARMA4  CL  test  sequence  with  5  dB 
of  additive  white  noise  using  LSMYWE  and  the  smallest  singular  value 
discarded.  The  MA  part  was  found  using  least  squares  identification. 
This  example  corresponds  to  the  example  in  Figure  (3.15a).  (a)  The  noisy 
sequence  and  (b)  the  estimated  model  impulse  response  and  the  original 
noiseless  sequence.  Model  order  (8,8). 
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IV.  LINEAR  MODELING  OF  ACOUSTIC 

TRANSIENTS 


A.  ACOUSTIC  DATA— GENERAL 

The  laboratory  generated  acoustic  transient  data  available  for  this  study  were 
generated  in  six  different  ways  using  ordinary  laboratory  objects.  The  data  records 
and  their  method  of  generation  are  summarized  in  Table  4.1.  Each  transient  will 
be  refered  to  by  the  object  used  or  the  action  performed  to  generate  that  particular 
transient  such  as  ‘book’  or  ‘slam’.  The  six  data  records  modeled  in  this  section  are 
shown  in  Figures  4.1a-f.  The  range  of  data  that  was  actually  modeled  is  listed  in 
Table  4.1.  The  sampling  rate  for  the  data  is  unknown  but  is  not  required  since  there 
is  no  need  to  infer  the  specific  characterics  of  the  acoustic  source  or  the  acoustic 
enviroment. 


Data  .Name 

Generation  Technique 

Indices  of  Modeled  Data 
(begin:end)  ■ 

Book 

Dropped  Book 

(55:4.54)  | 

Slam 

Slammed  Book 

(15:414) 

Plate 

Struck  Metal  Plate 

(5:704 )  i 

Ruler 

Book  Struck  with  Ruler 

(501:650) 

Wood 

Clapped  Wooden  Blocks 

(171:320)  1 

Wrench 

Dropped  Wrench 

(101:250)  ! 

TABLE  4.1:  Summary  of  Acoustic  Transient  Data  and  method  of  genera¬ 
tion. 


B.  ACOUSTIC  TRANSIENT  MODELING  RESULTS 

Iterative  prefiltering  and  correlation  domain  iterative  prefiltering  were  found  to 
yield  the  most  effective  time  domain  match  (in  terms  of  squared  error)  between  the 
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original  signal  and  the  impulse  response  of  the  estimated  model.  LSMYWE  with 
the  smallest  singular  value  set  to  zero,  as  in  Figures  3.13a,b,  was  used  to  find  the 
initializing  model  poles  for  the  iterative  algorithms.  Removing  the  smallest  singular 
value  cillowed  the  LSMYWE  algorithm  to  place  zeros  much  closer  to  the  unit  circle 
than  was  possible  when  all  singular  values  were  used.  Prony’s  method  and  the  Akaike 
MLE  algorithm  were  also  used  to  model  laboratory  tramsients.  However,  their  perfor¬ 
mance  was  generally  poor  and  they  will  not  appear  in  the  remainder  of  this  chapter. 
As  noted  in  Chapter  III,  the  problem  with  Prony’s  method  is  that  it  camnot  easily 
model  highly  resonant  frequencies,  that  is,  zeros  close  to  the  unit  circle,  in  the  pres¬ 
ence  of  even  small  amounts  of  noise  unless  a  large  number  of  excess  poles  are  used 
and  numerous  singular  values  are  discarded.  This  procedure  is  burdensome  when 
initializing  iterative  algorithms  which  do  not  require  these  excess  poles.  The  primary 
difficulty  with  Akaike  MLE  i.s  its  inability  to  handle  zeros  outside  the  unit  circle.  The 
initializing  model  zeros  were  found  using  Shank's  method  which  was  by  far  the  most 
robust  algorithm  for  finding  numerator  coefficients  of  the  methods  tested  in  Chapter 
Three.  In  addition  to  the  time  domain  plots  of  model  impulse  response,  a  pole-zero 
model  spectrum  and  pole-zero  plot  is  be  provided  for  each  transient  so  that  their 
differing  characteristics  can  be  observed.  All  spectra  were  generated  by  squaring  the 
FFT  magnitude  of  either  the  model  coefficients  or  the  signal  being  modeled.  Model 
order  was  chosen  based  on  the  best  educated  guess  of  the  author  in  accordance  with 
the  recommendations  in  Chapter  Ill.  and  augmented  by  trial  and  error. 

The  LSMYWE  model  used  to  initialize  the  iterative  prefiliering  algorithm  for 
the  Slam  transient  is  shown  in  Figures  4.2a  and  4.3a.  The  corresponding  iterative 
prefiltering  model  shown  in  Figures  4.2b  and  4.3b.  Figure  4.4a  shows  the  iterative 
prefiltering  model  spectrum.  Figure  4.4b  illustrates  a  characteristic  of  the  itera¬ 
tive  prefiltering  algorithm  that  was  observed  throughout  the  modeling  of  acoustic 
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transients.  Namely,  when  an  excess  number  of  model  parameters  are  used,  the  error 
between  the  model  impulse  response  and  the  signal  being  modeled  increases  dramat¬ 
ically,  mainly  at  the  beginning  of  the  signal.  The  subsequent  application  of  Shank’s 
method  is  effective  in  reducing  this  initial  error.  This  effect  is  shown  in  Figure  4.5a. 
In  general,  the  application  of  Shank’s  method  as  the  last  step  in  the  modeling  process 
was  found  to  reduce  the  the  mean  squared  error  between  the  original  signal  and  the 
model  impulse  response  to  some  degree  for  all  modeling  methods.  The  sensitivity  to 
excess  parameters  shown  by  iterative  prefiltering  does  not  affect  correlation  domain 
iterative  prefiltering.  Note  that  excess  model  zeros  cause  no  difficulties  in  Figure 
4.5b.  For  the  Book  transient,  two  modeling  trials  are  shown.  Figures  4.6  and  4.7 
show  the  best  time  domain  match  obtained  using  iterative  prefiltering  and  also  an 
LSM’t’WE.  Shank's  method,  respectively.  Although  the  LSNH  WE  model  does  not 
achieve  as  effective  an  impulse  response  match  as  iterative  prefiltering,  with  S\T)  it 
is  more  sensitive  to  the  low  energy',  high  frequency  component  present  in  the  Book 
transient  at  approximately  j.  The  best  model  of  Ruler  is  shown  in  Figure  4.S.  The 
model  spectra  for  Book  and  Ruler  appear  in  Figure  4.9. 

The  assumption  that  the  signal  being  modeled  is  a  system  impulse  response 
is  problematic  for  the  Plate,  Wood,  and  Wrench  signals  since  they  do  not  exhibit 
the  rapid  decay  usually  associated  with  an  impulse  response.  However,  it  is  still 
possible  to  achieve  a  resonable  time  domain  match  over  small  segments  of  each  signal. 
This  result  is  illustrated  in  the  models  of  in  Figures  4.10,  4.11.  4.13,  and  4.14.  The 
model  spectra  for  Wood  and  Wrench  arc  shown  in  Figure  4.15.  The  reason  that 
correlation  domain  iterative  prefiltering  was  used  for  the  Plate.  Wood,  and  Wrench 
signals  instead  of  standard  iterative  prefiltering  can  be  illustrated  by  comparing  the 
two  techniques  on  a  segment  of  the  Plate  signal.  Although  the  impulse  response  error 
of  the  two  models  in  Figures  4.10  and  4.11  are  nearly  identical.  Figure  4.12  shows 
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that  correlation  domain  iterative  prefiltering  clearly  outperformed  standard  iterative 
prefiltering  in  reproducing  the  spectrum  ^f  the  original  sequence.  In  f2w:t,  at  the  more 
suitable  model  order  of  (16,16),  iterative  prefiltering  would  not  converge  but  instead 
oscillated  in  a  region  of  convergence  for  any  model  order  over  (12,12).  Figure  4.16b 
shows  the  model  impulse  response  obtained  when  the  large  segment  of  Plate  shown 
in  Figure  4.16a  is  modeled  using  iterative  prefiltering.  Although  the  time  domain  and 
spectral  properties  of  the  model  relative  to  the  original  signal  are  considerable  poorer 
than  those  obtained  for  a  short  segment,  the  model  does  clearly  share  many  of  the 
features  of  the  original  signal. 

C.  ACOUSTIC  SIGNAL  MODELING  SUMMARY 

The  modeling  results  obtained  in  the  previous  section  of  this  Chapter  indicate 
the  possible  utility  of  pole-zero  modeling  algorithms  with  regard  to  modeling  tran¬ 
sient  signals.  Signals  with  decaying  narrowband  components  (e.g.  Slam.  Book,  and 
Ruler)  and  signals  with  substained  narrow'band  components  (e  g.  Plate.  Wood,  and 
Wrench  )  can  be  modeled  as  the  impulse  response  of  a  rational  linear  system.  Robust 
modeling  algorithms  are  available  which  can  effectively  deal  with  the  many  uncertain¬ 
ties  associated  with  real  world  signals.  .Although  the  goal  of  achieving  an  exact  time 
domain  match  between  the  original  signal  and  the  pole-zero  model  impulse  response 
was  not  realized  for  any  of  the  acoustic  signals  in  this  chapter,  in  all  cases  the  degree 
of  match  obtained  clearly  indicates  that  many  signal  characteristics  are  described  bv 
the  model. 
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Figure  4.1:  Laboratory  generated  acoustic  transient  data:  (a)  Slam  and 
(b)  Book. 
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Figure  4.1:  continued  Laboratory  generated  acoustic  transient  data:  (c) 
Ruler  and  (d)  Plate. 
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Figure  4.1:  continued  Laborotory  generated  acoustic  transient  data:  (e) 
Wood  and  (f)  Wrench. 
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Figure  4.2:  Modeling  the  transient  Slam — model  impulse  response  versus 
the  original  signal.  Model  order  (6,8).  (a)  LSMYi\'E  and  Shank's  method 
and  (b)  iterative  prefiltering  initialized  with  (a). 
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Figure  4.3:  Modeling  the  transient  Slam — model  pole-zero  plots.  Model 
order  (6,8).  (a)  LSMYWE  with  SVD  and  Shank's  method  and  (b)  iterative 
prefiltering  initialized  with  (a). 
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Figure  4.4:  Modeling  the  transient  Slam — model  spectrum  and  and  over¬ 
parametrized  iterative  prefiltering  model  impulse  response  versus  the  orig¬ 
inal  signal,  (a)  Model  spectrum  for  iterative  preiiltering  with  model  order 
(6,8)  and  (b)  iterative  prefiltering  for  model  order(6,12).  Note  how  excess 
model  parameters  cause  error  at  the  beginning  of  the  iterative  prefiltering 
model. 
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Figure  4.5.  N4odeiing  the  transient  Slun — model  impulse  response  versus 
the  original  signal,  overparameterized  modeling  effects,  (a)  The  applica¬ 
tion  of  Shanks  method  to  the  iterative  prefiltering  model  of  order  (6,12) 
reduces  the  error  at  the  beginning  of  tb*  model  impulse  response  and 
(b)  correlation  domain  iterative  prefiltcring  for  model  order  (6,12)  can 
effectively  use  the  excess  model  parameterf^. 
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Figure  4.6:  Modeling  the  transient  Book — model  impulse  response  versus 
the  original  signal.  Model  order  (6,6).  (a)  LSMYWE  with  SVD  and 

Shank’s  method  and  (b)  iterative  prefiltering  initialized  with  (a).  Note 
the  sensitivity  of  LSMYWE  with  SVD  to  the  high  frequency  components 
present. 
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Figure  4.7:  Modeling  the  transient  Book — model  pole-zero  plots.  Model 
order  (6,6).  (a)  LSMYWE  with  SVD  and  Shank's  method  and  (b)  iterative 
prefiltering  initialized  with  (a).  Note  the  sensitivity  of  LSMYWE  with 
SVD  to  the  high  frequency  pole  pair. 


77 


-  original  signal 

-  model  impulse  response 


X  model  poles 
o  model  zeros 


-1 . 5  ^ - 

z-plane 


(b) 


Figure  4.8:  Modeling  the  transient  Ruler — Model  order  (6,16).  (a)  Iter¬ 
ative  prefiltering  model  impulse  response  versus  the  original  signal  and 
(b)  the  corresponding  model  pole-scro  plot.  Note  that  the  two  beating 
sinusoids  have  been  effectively  modeled. 
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Figure  4.9:  Modeling  the  transient  Book  and  Ruler — model  ipectrums. 
(a)  The  iterative  prefiltering  model  of  order  (6,6)  for  the  transient  Book 
and  (b)  the  iterative  prefiltering  model  of  order(6,12)  for  the  transient 
Ruler. 


79 


-  original  signal 

-  model  impulse  response 


0  50  100  150 


n 

X  model  poles 
o  model  zeros 


-2  ^ - - 

z-plane 

Figure  4.10:  Modeling  the  transient  Plate — Model  order  (12,12).  (a)  It¬ 
erative  prefiltering  model  impulse  response  versus  the  original  signal  and 
(b)  the  corresponding  model  pole-xero  plot. 
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Figure  4.11;  Modeling  the  transient  Plate — Model  order  (16,16).  (a)  Cor¬ 
relation  domain  iterative  prefiltering  model  impulse  response  versus  the 
original  signal  and  (b)  the  corresponding  model  pole-zero  plot. 
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Figure  4.12;  Modeling  the  transient  Plate — ^model  ipectrums  versus  the 
original  signal  spectnims.  (a)  The  iterative  prefiltering  model  of  order 
(12,12)  and  (b)  the  correlation  domain  iterative  prefiltering  model  of  order 
(16,16) 
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Figure  4.13:  Modeling  the  transient  Wood — Model  order  (16,16).  (a) 

Correlation  domain  iterative  prefiltering  model  impulse  response  versus 
the  original  signsd  and  (b)  the  corresponding  model  pole-zero  plot. 
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Figure  4.14:  Modeling  the  transient  Wrench — Model  order  (8,8).  (a)  Cor¬ 
relation  domain  iterative  prefiltering  model  impulse  response  versus  the 
original  signal  and  (b)  the  corresponding  model  pole-zero  plot. 
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Figure  4.15:  Modeling  the  transient  Wood  and  Wrench — model  spectrums. 
(a)  The  iterative  prefiltering  model  of  order  (16,16)  for  the  transient  Wood 
and  (b)  the  iterative  prefiltering  mod*l  of  order(8,6)  for  the  transient 
Wrench. 


85 


1000 


ufoTism*®’-  ‘l>«  ‘r.n.ient  Pl.t.  ov.r  the  Urge  eegment 

(450^250)  using  iterative  prefiltering,  model  order  (16,16).  (a)  The  oriir- 
inal  Plate  segment  and  (b)  the  model  impulse  response.  * 


86 


V.  CONCLUSIONS 


A.  PERFORMANCE  COMPARISON  SUMMARY 

The  modeling  of  signals  as  the  impulse  response  of  a  linear  pole-zero  system  is 
an  important  tool  in  digital  signal  processing.  One  natural  area  for  applying  such 
an  approach  is  in  the  modeling  of  tr£insient,  impulse  response-like  waveforms.  The 
specific  approach  taken  in  this  thesis  was  to  determine  which  pole- zero  modeling  al¬ 
gorithms  are  most  suited  to  modeling  complex,  real  world  transient  waveforms.  The 
modeling  criterion  emphasized  is  to  obtain  the  best  (least  squares  error)  time  domain 
match  between  the  model  impulse  response  and  the  original  signal.  This  criterion 
was  adopted  because  it  provides  a  degree  of  signal  characterization  that  is  a  step  be¬ 
yond  normal  power  spectrum  estimation.  Indeed,  the  strength  of  pole-zero  models  is 
their  ability  to  describe  not  only  the  resonances  present  (model  poles),  but  also  how 
these  resonances  are  related  (model  zeros).  Because  of  the  widely  varying  character¬ 
istics  anticipated  for  real  world  signals,  a  key  evaluation  criterion  for  any  modeling 
algorithm  is  robustness.  In  particular,  algorithms  must  be  effective  in  the  presence  of 
signal  degrading  effects  like  noise  and  model  degrading  effects  such  as  unknown  model 
order.  The  performance  of  several  selected  algorithms  were  compared  for  known  im¬ 
pulse  response  test  sequences  in  Chapter  III.  The  modeling  experience  gained  in  these 
experiments  was  then  applied  to  modeling  laboratory  generated  acoustic  transient 
datain  in  Chapter  IV'  as  a  test  of  ‘real  world'  effectiveness. 

Four  basic  algorithms  were  chosen  for  comparison;  Prony's  method,  the  least 
squares  modified  Yule-Walker  equations  (LSMYWE),  iterative  prefiltering,  and  the 
.Akaike  maximum  likelihood  estimator  (MLE).  An  algorithm  which  is  an  extension  of 
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iterative  prefiltering  into  the  correlation  domain  was  also  presented.  For  those  meth¬ 
ods  in  which  the  model  poles  and  zeros  are  determined  separately  (Prony’s  method 
and  LSMYWE),  four  methods  for  determining  model  zeros  (i.e.,  transfer  function  nu¬ 
merator  coefficients)  were  considered;  the  upper  partition  of  Prony's  method,  spectral 
factorization,  Durbin’s  method,  and  Shank’s  method.  The  major  conclusions  of  the 
algorithm  performance  comparisons  conducted  in  Chapter  III  and  Chapter  I\  are  as 
follows; 

1.  .Algorithms  that  are  unable  to  model  zeros  outside  the  unit  circle  (Durbin's 
method,  Akaike  NILE)  have  limited  versatility  when  modeling  arbitrary  tran¬ 
sient  waveforms.  All  the  acoustic  transients  in  Chapter  IV'  required  a  non¬ 
minimum  phase  model  to  obtain  the  best  time  domain  match, 

2.  The  most  robust  and  effective  method  for  finding  zeros  that  gives  the  best  least 
squares  time  domain  match  was  found  to  be  Shank's  method.  In  fact,  applving 
Shank's  method  as  a  last  step  improved  the  final  model  of  all  algorithms  to 
some  degree.  The  upper  partition  of  Prony's  method  and  spectral  factorization 
are  not  very  useful  because  of  their  extreme  sensitivity  to  noisy  or  time  shifted 
signals. 

3.  Prony  s  method  and  Akaike  MLE  have  difficulty  modeling  signals  in  which 
additive  noise  is  present.  Even  small  amounts  of  additive  noise  causes  a  dramatic 
loss  of  pole  frequency  resolution  for  Prony's  method  The  use  of  excess  poles 
and  singular  value  decomposition  were  found  to  be  effective  in  overcoming  these 
effects  but  these  methods  depend  on  a  knowledge  of  correct  model  order.  .Also, 
unlike  spectrum  estimation,  excess  poles  must  be  retained  for  time  domain 
matching.  This  is  in  direct  contrast  to  the  parsimonious  use  of  model  parameters 
normally  provided  by  a  pole-zero  model.  The  difficulties  encountered  with  the 
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Akaike  MLE  were  computational  in  nature;  for  noisy  signals  the  Akaike  MLE 
algorithm  usually  terminated  prematurely  when  either  a  non-nxinimum  phase 
model  was  encountered  during  an  iteration  or  when  a  numerical  overflow  was 
induced  by  an  unstable  model  estimate. 

4.  LSMV’WE  with  singular  value  decomposition  and  iterative  prefiltering  (includ¬ 
ing  the  correlation  domain  version)  were  found  to  be  the  most  effective  algo¬ 
rithms  for  modeling  a  signal  when  additive  noise  is  present.  If  the  true  model 
order  of  a  system  is  unknown,  it  is  best  to  discard  only  one  singular  value. 
.Almost  all  the  resolution  gain  occurs  with  the  first  singular  value.  Discarding 
additional  singular  values  is  intended  to  reduce  the  variance  of  excess  poles  but 
it  will  cause  poles  to  be  biased  if  all  model  poles  are  necessary  for  signal  model¬ 
ing.  Both  of  these  methods  demonstrated  the  consistent  ability  to  model  poles 
very  close  to  the  unit  circle.  This  capability  was  es.sential  when  modeling  the 
acoustic  transients  used  in  this  thesis. 

Combining  the  above  observations  leads  to  our  recommended  strategy  for  mod¬ 
eling  an  arbitrary  transient  signal.  First,  use  LSM’t'WE  with  one  singular  value  re¬ 
moved  and  Shank’s  method  to  find  an  initial  model  estimate.  Next,  use  this  estimate 
to  initialize  either  iterative  prefiltering  or  correlation  domain  iterative  prefiltering.  Fi¬ 
nally.  if  desired,  apply  Shank's  method  to  optimize  the  time  domain  fit  of  the  model 
zeros. 

B.  RECOMMENDATIONS  FOR  FUTURE  STUDY 

The  results  of  Chapter  IV  indicate  that  there  are  pole- zero  modeling  algorithms 
available  that  are  sufficiently  robust  to  be  useful  for  modeling  many  complex,  real 
world  transients.  A  number  of  issues  regarding  the  pole-zero  modeling  of  transient 
signals  and  the  application  of  such  models  require  further  study.  These  issues  include; 
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1.  The  degree  to  which  the  noise  performance  of  correlation  domain  iterative  pre¬ 
filtering  may  be  improved  by  introducing  noise  compensation  of  the  autocorre¬ 
lation  sequence  requires  study. 

2.  A  study  of  the  effectiveness  of  model  order  determination  techniques  when  ap¬ 
plied  to  transient  modeling  would  facilitate  more  effective  use  of  transient  mod¬ 
eling  techniques. 

3.  The  effectiveness  of  iterative  prefiltering  and  correlation  domain  prefiltering  a.s 
a  spectral  estimation  technique  should  be  explored  further. 

4.  The  ability  of  pole-zero  models  to  describe  the  time  domain  characteristics  of 
a  signal  could  aid  in  the  detection  of  signals.  Such  an  application  needs  to  be 
persued. 

5.  A  study  of  the  structural  relationship  of  pole- zero  models  to  the  specific  sy.stpms 
that  generate  transients  may  increase  the  usefulness  of  such  models. 


90 


T 


f 


REFERENCES 

1.  George  E.  P.  Box  and  Gwilyn  M.  Jenkins,  Time  series  analysis,  forecasting  and 
control,  rev.  ed.,  Holden  Day,  San  Francisco,  1976. 

2.  John  M.  Gottman,  Time-series  analysis,  a  comprehensive  introduction  for  social 
scientists,  Cambridge  University  Press,  1981. 

3.  Mike  West  and  Jeff  Harrison,  Bayesian  forecasting  and  dynamic  models.  Springer 
Verlag,  New  York,  1989. 

4.  Steven  M.  Kay,  Modem  spectral  estimation:  theory  and  applications.  Prentice 
Hall  Inc.,  Englewood  Cliffs,  New  Jersey,  1988. 

5.  James  A.  Cadzow,  “Spectral  estimation:  an  overdetermined  rational  model  equa¬ 
tion  approach,”  Proc.  IEEE,  Vol.  70,  No.  9,  September  1982. 

6.  S.B.  Kesler,  ed.,  ’’Modern  spectrum  analysis  11,”  selected  reprint  series  of  IEEE 
Acousl.  Speech  and  Sig.  Proc.  Society,  IEEE  Press,  1986. 

7.  N.K.  Sinha  and  B.  Kuszta,  ”Modeling  and  identification  of  dynamic  systems.” 
Van  Nostrand  Reinhold  Company,  New  York.  1983. 

8.  Pieter  Eykoff,  “System  identificationf'  John  Wiley  and  Sons,  New  York.  1974. 

9.  Lennart  Ljung,  System  identification:  theory  for  the  user.  Prentice-Hall  Inc.. 
Englewood  Cliffs,  New  Jersey,  1987. 

10.  Torsten  Soderstrum  and  Petre  Stoica,  System  identification.  Prentice  Hall.  Lon¬ 
don,  1989. 

11.  Lawrence  R.  Rabiner  and  Ronald  W’. Schafer,  Digital  signal  processing  of  speech 
signals.  Prentice-Hall  Inc.,  Englewood  Cliffs,  New  Jersey.  1978. 

12.  R.  Genesio  and  M.  Milanese,  “A  note  on  the  derivation  and  use  of  reduced-order 
models,"  IEEE  Trans.  Automat.  Contr.,  V’ol.  AC-21,  pp.  118-122,  February  1976. 

13.  Cyril  M.  Harris,  ed..  Shock  and  vibration  handbook.  McGraw-Hill  Book  Companv, 
New  York,  1988. 

14.  D.  Brown,  G.  Carbon,  and  K.  Ramsey,  "Survey  of  excitation  techniques  applicable 
to  the  testing  of  automotive  structures,”  Soc.  Automotive  Engineers,  SAE  Paper 
77029.  1977. 

15.  C.E.  Baum,  "Emerging  technolo©'  for  transient  and  broad  band  analysis  and 
synthesis  of  antenna  scatterers,"  Proc.  IEEE,  VoX.  64,  pp.  1598-1616,  1976. 

16.  Dale  G.  Stone,  “Wavelet  estimation,"  Proceed.  IEEE,  Vol.  bf  72,  No.  10.  pp. 
1394-1402.  October  1984. 

17.  Micheal  A.  Van  Blaricum  and  Raj  Mittra,  "Problems  and  solutions  associated 
with  Prony’s  method  for  processing  transient  data,”  IEEE  Trans,  on  An/,  and 
Prop.,  Vol.  AP-28,  No.  1,  pp.  174-182,  January  1978. 

18.  T.L.  Henderson.  "Geometric  methods  for  determining  system  poles  from  transient 
response,”  IEEE  Trans,  on  Acoust.,  Speech,  and  Sig.  Proc.,  Vol.  ASSP-29,  No.  5. 
pp.  982-988,  October  1881. 


91 


I 


19.  Sophocles  J.  Orfanides,  "Pole  retrieval  by  eigenvector  methods,”  IEEE  Inti  Conf. 
Acoust.,  Speech,  and  Sig.  Proc.,  pp.  1505-1508,  1987. 

20.  D.L.  Brown,  et.al.,  “Parameter  estimation  techniques  for  modal  analysis,”  Soc. 
Automotive  Engineers,  SAE  Paper  790221,  1979. 

21.  Ramdas  Kum^u•esan  and  Donald  W.  Tufts,  ’’Estimating  the  pau-ameters  of  ex¬ 
ponentially  damped  sinusoids  and  pole-zero  modeling  in  noise.”  IEEE  Trans. 
Acoust.  Speech  Sig.  Proc.,  Vol.  ASSP-30,  pp.  833-840,  December  1982. 

22.  James  H.  McClellan,  “Parametric  signal  modeling,”  in  Advanced  topics  in  signal 
processing,  Jae  S.  Lim  and  Alan  V.  Oppenheim,  eds.,  pp.  1-57,  Prentice  Hall. 
Englewood  Cliffs,  New  Jersey,  1988. 

23.  Murali  Tummula,  “Iterative  method  for  ARM  A  least  squares  identification,”  sub¬ 
mitted  for  publication. 

24.  Steiglitz,  K.  and  McBride,  L.E.,  “A  technique  for  the  identification  of  linear 
systems,”  IEEE  Trans.  Automat.  Contr.,  Vol.  AC-10,  pp.  461-464.  October  1965. 

25.  Alan  G.  Evans  and  Robert  Fischl,  “Optimal  least  squares  time-domain  synthesis 
of  recursive  digital  filters,”  IEEE  Trans.  Audio  Elect oacoust.,  Vol.  AU-21.  No.  I. 
pp.  61-65,  February  1973. 

26.  Jackson,  L.B.,  “Digital  Filters  and  SignsJ  Processing.”  Second  Edition.  Kluwer 
Academic  Publishers.  Boston,  1989. 

27.  Ira  S.  Konvalinka  and  Miroslav  R.  Matausek.  “Simultaneous  estimation  of  poles 
and  zeros  in  speech  analysis  and  ITIF — itrative  inverse  filtering  algorithm."  IEEE 
Trans.  Acoust.,  Speech  and  Signal  Process.,  V'ol.  ASSP-27,  No.  5.  pp.  485-492. 
October  1979. 

28.  M.  Isabel  Ribeiro  and  Jose  M.F.  Moura,  “Dual  algoritm  for  .AR.M.A  spectrum 
estimation.”  IEEE  Conf.  Acoust.,  Speech  and  Signal  Process.,  pp.  2081-2084. 

29.  Jitendra  K  Tugnait.  “Identification  of  linear  stochastic  systems  via  .secuiid-  and 
fourth-order  cumulant  matching.”  IEEE  Trans.  Inform.  Theory.  \'oI.  lT-33.  No 
3.  pp.  39.3-407,  .May  1987. 

30.  Georgios  B.  Giannakis  and  Jerry’  M  Mendel,  “Stochastic  realization  of  non¬ 
minimum  phase  systems,”  Proc.  American  Control  Conf  (ACC),  pp.  1254  1259. 
Seattle,  WA,  June  1986. 

.31 .  H.  Akaike.  “Maximum  likelihood  identification  of  gaussian  autoregressive  moving- 
average  models.”  fliometnca,  Vol.  60.  pp.  255-265.  .August  1973 

.32.  Craig  F.  Ansley,  “An  algorithm  for  the  exact  likelihood  of  a  mixed  autoregressive 
moving  average  process,”  Biometrika,  Vol.  66,  No.  1.  pp.  .59-65.  1979. 

33.  N.H.  Judell,  “Maximum  likelihood  parameter  estimation  for  signals  with  rational 
spectra.”  M.S.  Thesis,  Univ.  of  Rhode  Island,  1983. 

34  Leland  B.  Jackson,  Jianguo  Huang,  and  Kevin  P.  Richards,  “.AR.  ARMA.  and 
AR-in-noise  modeling  by  fitting  windowed  correlation  data.”  IEEE  Conf  .Acoust.. 
Speech  and  Signal  Process.,  pp.  2039-2042,  1987. 

.35  R.  Moses  et.  al.,  “An  efficient  linear  method  of  ARMA  spectral  estimation."  lEF.E 
Conf.  Acoust.,  Speech  and  Signal  Process.,  pp.  2077-2080,  1987. 


92 


36.  Donald  W.  Tufts  and  Ramdas  Kumarestin,  “Estimation  of  frequencies  of  multi¬ 
ple  sinusoids:  making  linear  prediction  perform  like  m^lximum  likelihood,”  Proc. 
IEEE,  Vol.  70,  No.  9,  pp.  975-989,  September  1982. 

37.  Gerard  Alengrin  and  Josiane  Zerubia,  “A  method  to  estimate  the  parameters 
of  an  ARMA  model,”  IEEE  Trans.  Automat.  Contr.,  Vol.  AC-32,  No.  12,  pp. 
1113-1115,  December  1987. 

38.  Haykin,  S.,  “Modern  Filters,”  Macmillan  Publishing  Company,  New  York.  1989. 

39.  James  A.  Cadzow,  “Spectral  analysis,”  in  D.  F.  Elliot,  ed..  Handbook  of  digital 
signal  processing,  .Academic  Press,  Inc  ,  1987. 

40.  Therrien,  C.W.,  “Discrete  random  signals  and  statistical  signal  processing,"  Un¬ 
published  Manuscript,  2  April  1990. 

41.  J.  Durbin,  “Efficient  estimation  of  parameters  in  moving- average  models." 
Biometrika,  Vol.  46,  pp.  306-316,  1959. 

42.  John  L.  Shanks,  “Recursion  filters  for  digital  processing,"  Geophysics,  \'ol.  32. 
No.  1,  pp.  33-51,  February  1967. 

43.  Steiglitz,  K.,  “On  the  simultaneous  estimation  of  poles  and  zeros  in  speech  anal¬ 
ysis,”  IEEE  Trans.  Acoust.  Speech  Signal  Proc.,  Vol.  .ASSP-25.  pp.  229-234.  June 
1977. 

44.  Torsten  Soderstrom,  “Identification  of  stochastic  linear  systems  in  the  presence 
of  input  noise,"  Automatica,  Vol.  17,  No.  5.  pp.  713-72.5.  1981. 

45.  Sir  Ronald  .A.  Fisher.  “Statistical  methods  and  scientific  inference."  Third  Ed.. 
Hafner  Press.  .A  Division  of  Macmillam  Publishing  Co..  Inc..  New  Yoik.  1973. 

46.  K.J.  Astrom.  P.  Hagander.  and  J  Sternbv,  “Zeros  of  sampled  systems."  .Automat- 
ica.  \ol.  20.  pp.  31-38.  1984. 

47.  J..\l.  .Mendel.  “Optimal  Siesmtc  Deconvolution."  Academic  Press.  .New  York. 
1983. 

4^.  .Albert  Benveniste.  .Maurice  Goursat,  and  Gabriel  Ruget.  “Robust  identification 
of  a  nonminimum  phase  system:  blind  adjustment  of  a  linear  equalizer  in  data 
communications."  IEEE  Trans.  Automat.  Contr..  \bl.  .AC-25.  No.  3,  pp.  3S5  39.S. 
June  1980. 

49  J.S.  Freudenburg  and  D.P.  Looze.  “Sensitivity  reduction,  nonminimum  phase 
zeros,  and  design  tradeoffs  in  single  loop  feedback  systems,"  IEEE  Inti.  Conf. 
Dec.  and  Contr.,  pp.  625-630.  December  1983. 

50.  John  J.  Kormylo  and  Jerry  M  Mendel.  “Identifiability  of  nonminimum  phase 
linear  stochastic  systems,”  IEEE  Trans.  Automat.  Control.  Vol.  .A('-28.  .No,  12. 
pp.  1081-1089.  December  1983. 

51.  Vijay  K.  Arya  and  H.D.  Holden,  “Deconvolution  of  seismic  data — an  overview." 
IEEE  Trans.  Geosci.  Electron.,  Vol.  GEi-16,  No.  2.  pp.  95-98.  .April  1978. 

52.  Hirotugo  Akaike.  “.A  new  look  at  statistical  model  identification,"  IEEE  Tran.<t. 
Automat.  Contr.,  Vol.  AC-19,  No.  6.  pp.  716-723.  December  1974. 

53  E.J.  Hannan,  “The  estimation  of  the  order  of  an  ARMA  process,”  .Annah  of  Stat.. 
Vol.  8.  .No.  5.  pp.  1071-1081.  1980. 


93 


54.  Rangasami  Kashyap,  “Optimal  choice  of  AR  and  MA  parts  in  autoregressive 
moving  average  models,”  IEEE  Trans.  Patt.  Analy.  and  Mach.  Intelligence.  Vol. 
PAMI-4,  No.  2,  pp.  99-104,  March  1982. 

55.  Steven  M.Kay,  “The  effects  of  noise  on  the  autoregressive  spectral  estimator.” 
IEEE  Trans,  on  Acoust,  Speech,  and  Signal  Proc.,  Vol.  ASSP-27,  No.  5,  pp.  47S- 
485,  October  1979. 

56.  Leland  B.  Jackson,  et.  al.,  “Frequency  estimation  by  linear  prediction."  Proc. 
Inti.  Conf.  Acoust.,  Speech,  Signal  Proc.,  pp.  352-356,  1978. 

57.  Steven  M.  Kay,  “Noise  compensation  for  autoregressive  spectral  estimates."  IEEE 
Trans.  Acoust.,  Speech,  Signal  Process..  Vol.  ASSP-28.  No.  3,  pp.  292-303.  June 
1980. 

58.  Luis  Vergara- Dominguez,  “New  insights  into  the  high-order  Yule- Walker  equa¬ 
tions,”  IEEE  Trans.  Acoust.,  Speech.  Signal  Process..  Vol.  .\SSP  38.  No.  9.  pp. 
1649-1651,  September  1990. 

59.  James  .\.  Cadzow  and  Behshad  Baseghi.  “Data  adaptive  .4R.\1.\  modeling  of  tirn<- 
series,"  Proc.  Inti.  Conf.  Acoust.,  Speech.  Signal  Process.,  pp.  256-261.  1982 

60.  Sriniv£isan,  T.,  Swanson.  D.C..  Symons,  F.W..  “ARM.^  model  order/data  length 
tradeoff  for  specified  frequency  resolution,"  Proc.  IC.\SSP.  198.}.  Paper  38.2. 

61 .  Petre  Stoica  and  Torsten  Soderstrom.  “The  Steiglitz-.McBridp  identification  alco- 
rithm  revisited — convergence  analysis  and  accuracy  aspects."  IEEE  Tran.^  .lu- 
tomat.  Contr..  Vol.  AC-26,  No.  3,  pp.  712-717.  June  1981. 

62.  Donald  W.  Tufts  and  R.  Kumaresan.  “Improved  spectral  resolution."  Proc.  IEEE. 
Vol.  68.  No.  3.  pp.  419-420.  March  1980. 


94 


